!!****************************************************************************
!! FUNCALLER: i.  solves for default model
!!            ii. generates simulations and computes MC statistics
!!****************************************************************************

REAL(8) FUNCTION func(x0)
    
USE GLOBAL
USE OMP_LIB

implicit none

REAL(8),INTENT(IN)::x0(6)

REAL(8)::v2(na,nb2),c1(na,nb2),uc1(na,nb2)
REAL(8)::wmax,p0max,h0max,gn0max,ct0max,m0max,tau0max,ww0max,f0max,b_next_global
REAL(8)::value1,value2,objective_u_fe,objective_u_unemp,objective_u_tax
INTEGER  i,i_b, i_a, i_default_global,now_time
INTEGER::index_b,new_index_b,number_of_threads
INTEGER(4)::hours2,mins2,index_opt 
REAL(8)::a_initial,b_initial,udef,u_value,indicator_external
REAL(8)::stats(6),xlb(6),xub(6),penalty,time,ct_matrix_new(nb, na)
EXTERNAL udef, objective_u_fe,objective_u_unemp,objective_u_tax
REAL(8)::end_time_loop,start_time_loop      
REAL(8)::old_value,tau_value2,new_value,tau_value 

iterfun = iterfun+1

beta =x0(1)
chi  =x0(2)
icost=x0(3)
psi1 =x0(4)
psi2 =x0(5)
wbar =x0(6)


indicator_external =1 !=0 compute initial guess, =1 reads initial guess from external files, >5 reads value functions and prices for finer grids
if (indicator_external < 0.5) then
        call compute_bgrid
        call INITIAL       
      
        open (10, FILE='graphs\v.txt',STATUS='replace')
        open (11, FILE='graphs\default.txt',STATUS='replace')
        open (12, FILE='graphs\x.txt',STATUS='replace')
        open (13, FILE='graphs\gn.txt',STATUS='replace')
        open (14, FILE='graphs\b_next.txt',STATUS='replace')
        open (16, FILE='graphs\ctt.txt',STATUS='replace')

        do i_b = 1,nb
            do i_a = 1,na 
 
                a_initial = agrid(i_a)
                b_initial = bgrid(i_b)

                i_default_global = 2 !defaults                 
                
                old_value=-10d33

                DO i=1,ntau
                    tau_value2 = taugrid(i)

                   new_value = objective_u_tax(0d0,a_initial,b_initial,i_default_global,tau_value2)
                   if (new_value >= old_value) then
                      index_opt = i
                      old_value = new_value
                  end if
                end do

                tau_value = taugrid(index_opt)
                tau_global = tau_value

                u_value = objective_u_tax(0d0,a_initial,b_initial,i_default_global,tau_value)                 
                vaut_matrix(i_b, i_a) = u_value-udef(a_initial)                    

                i_default_global = 1 !repays                 
                !b_next_global=(1d0-lambda)*b_initial
                b_next_global=bmin 

                IF (welfare_case.EQ.1) b_next_global=(1d0-lambda)*b_initial
                    
                old_value=-10d33

                DO i=1,ntau
                    tau_value2 = taugrid(i)

                   new_value = objective_u_tax(b_next_global,a_initial,b_initial,i_default_global,tau_value2)
                   if (new_value >= old_value) then
                      index_opt = i
                      old_value = new_value
                  end if
                end do

                tau_value = taugrid(index_opt)
                tau_global = tau_value

                u_value = objective_u_tax(b_next_global,a_initial,b_initial,i_default_global,tau_value)                 

                x_matrix_new(i_b, i_a)=x_global
                gn_matrix_new(i_b, i_a)=gn_global
                ct_matrix_new(i_b, i_a)=ctrad_global
                v0_matrix(i_b, i_a) = u_value          

                v_matrix(i_b, i_a) = MAX(vaut_matrix(i_b, i_a), v0_matrix(i_b,i_a))

                 if (vaut_matrix(i_b, i_a) <= v0_matrix(i_b, i_a)) then
                   default_decision(i_b, i_a) = 1
                    b_next_matrix(i_b, i_a) = b_next_global
                 else
                   default_decision(i_b, i_a) = 2
                    b_next_matrix(i_b, i_a) = 0d0
                 end if              
                 
                 WRITE(10, '(F18.10, X, F18.10, X, F15.10)') v_matrix(i_b, i_a),MAX(-100d0,v0_matrix(i_b, i_a)),MAX(-100d0,vaut_matrix(i_b,i_a))
                 WRITE(11, '(F18.10)') defaultgrid(default_decision(i_b, i_a))
                 WRITE(12, '(F18.10)') x_matrix_new(i_b, i_a)
                 WRITE(13, '(F18.10)') gn_matrix_new(i_b, i_a)
                 WRITE(14, '(F15.11)') b_next_matrix(i_b, i_a)
                 WRITE(16, '(F18.10)') ct_matrix_new(i_b, i_a)
                 
                 IF (psi1.LT.500d0) THEN        !default vs. risk-free debt economy         
                     q_matrix(i_b, i_a) = 0d0                 
                     q_matrix_nodef(i_b, i_a) = 0d0 
                 ELSE
                     q_matrix(i_b, i_a) = payoff/ (lambda + r)
                     q_matrix_nodef(i_b, i_a) = payoff/ (lambda + r)
                 ENDIF
                 
            end do
      end do
  
         CLOSE(10)
         CLOSE(11)
         CLOSE(12)
         CLOSE(13)
         CLOSE(14)
         CLOSE(16)
         
        PRINT*, 'No Initial Guess Uploaded'
        PRINT*,' '
elseif (indicator_external > 5) THEN    !reads value functions and prices computed for finer grids
    call compute_bgrid
    CALL INITIAL_V
else                                    !reads initial guess from external files

open (10, FILE='graphs\v.txt')
open (16, FILE='graphs\q_paid.txt')
open (17, FILE='graphs\b_next.txt')
open (23, FILE='graphs\bounds.txt')

      READ(23, '(F15.11, X, F15.11)') bmin, bmax
      
      call compute_bgrid

      do i_b = 1,nb
      do i_a = 1,na

             READ(10, '(F18.10, X, F18.10, X, F18.10)') v_matrix(i_b, i_a), &
                     v0_matrix(i_b, i_a), vaut_matrix(i_b, i_a)
             READ(16, '(F15.11, X, F15.11)') q_matrix(i_b, i_a), q_matrix_nodef(i_b, i_a)
             READ(17, '(F15.11)') b_next_matrix(i_b, i_a)
             if (vaut_matrix(i_b, i_a) > v0_matrix(i_b, i_a)) then
                 default_decision(i_b, i_a) =2
             else
                 default_decision(i_b, i_a) =1
             end if

        end do
        end do
CLOSE(10)
CLOSE(16)
CLOSE(17)
CLOSE(23)

PRINT*, 'Saved Initial Guess Uploaded'

end if 

!solves for baseline model
call iterate

!computes policy functions, price schedule and value functions using finer grid
call pol_fcns

!uploads policy functions, price schedule and value functions using finer grid
!CALL INITIAL_POL

!generates simulations and compute MC statistics
call MC_SIM(stats)

func = ((-stats(1)-22.8d0)/22.8d0)**2d0 +((stats(2)-18.1d0)/18.1d0)**2d0 +((stats(3)-1.1d0)/1.1d0)**2d0 +((stats(4)-1.4d0)/1.4d0)**2d0+((stats(5)-1.11d0)/1.11d0)**2d0+((stats(6)-0.9815d0)/0.9815d0)**2d0

call system_clock(count=now_time)      !Option B
time  = (now_time-start_time)/rate_time/MAX(1,iterfun-1)
hours2 = INT(FLOOR(time/3600))
mins2  = FLOOR((time/3600-(hours2))*60d0)
PRINT*,'Iter = ',iterfun, 'Time per iteration',hours2,'hs',mins2,'mins'


PRINT*, 'MODEL SOLVED AND MC STATS COMPUTED'
RETURN

END FUNCTION func
!******************************************************************************************
REAL(8) FUNCTION objective_u_fe(x,a_initial2,b_initial2,i_default_global2,tau_value)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,a_initial2,b_initial2,tau_value
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::bpval,wval,tauval,ct_global
REAL(8)::p0,m0,h0,q,b
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval
REAL(8)::aux_lhs,balt,aux0alt,minc,maxc,maxc2,maxc3
REAL(8)::term1,taxbar,bound1,bound2,new_value,old_value
REAL(8)::old_valuemax,sign_new,sign_old,sign_two,ctemp2,difc

REAL(8)::objective_u_feold,gn_global1old,x_global1old,ct_global1old 
INTEGER(4)::i_num_star,i_num_star2,index_optmax,i_num
REAL(8)::q2,aux02,p02,gnval2 ,temp1,temp2
INTEGER,PARAMETER::num=100
REAL(8)::ctemp(num),tolval

INTEGER(4)::i,jj,im,iax,feas,error,niter2,ftype,try,index_opt
LOGICAL::flag00
EXTERNAL q_fun

objective_u_feold = -1d10
i_num_star=0
tolval = 0d0

b0 = b_initial2
a0 = DEXP(a_initial2)

b = x
q = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0           !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))

ctval= a0+aux0
ct_global = ctval   

IF (ctval.LE.0d0) THEN
    objective_u_fe= huge_n +1d3*(0d0-ctval)            
    if (choice_glob.EQ.1) THEN
        gn_global1 = 0d0
        x_global1  = 0d0
        ct_global1 = 0d0
    ENDIF
    RETURN
ENDIF

!case 1: frictionless: no binding tax cap
!new with lumpsum taxes
IF (taxrate.GT.1d0) THEN    !lump-sum taxes
    try   = 100
ELSE                        !income tax rate
    try = 1
ENDIF

IF (chi.EQ.0d0) THEN
    cnval=1d0
    gnval = 0d0
ELSE
    term1 = (1d0-chi)/chi*(1d0-omegac)
    term1 = term1**(1d0/sigma)
    cnval = term1/(1d0+term1)
    gnval = 1d0-cnval
ENDIF

ctval = a0+aux0
hval  = 1d0


p0   = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
wval = alpha*p0
tauval = p0*gnval-aux0

2046 taxbar = tau_value*(a0+in*p0)

2047 IF (taxrate.LE.1d0) THEN     
    flag00 = ((p0.LE.0d0).OR.(wval.LT.wbar).OR.(gnval.GT.1d0).OR.(ctval.LE.0d0).OR.(tauval.GT.taxbar)) 
ELSE
    flag00 = ((p0.LE.0d0).OR.(wval.LT.wbar).OR.(gnval.GT.1d0).OR.(ctval.LE.0d0)) 
ENDIF

IF (.NOT. flag00) THEN
    pres  = p0
    hres  = hval
    ctres = ctval
    taures= tauval
    wwres = wval
    gnres = gnval
    feas  = 1d0
ELSE
    ctres   = 0d0    
    feas    = 0d0
    taures  = tauval
    hres    = 0d0
    pres    = p0
    gnres   = gnval
    wwres   = wval
ENDIF

!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(wwres.LT.wbar)*(wbar-wwres)+(ctres.LE.0d0)*(0d0-ctres))
!without lumpsum taxes
IF (tau_value.LE.1d0) penalty = penalty +1d6*(taures.GT.taxbar)*(taures-taxbar)
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GT.1d0)*(gnres-1d0)+(gnres.LT.0d0)*(0d0-gnres))

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres  = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres  = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

cres = cres - (cnres.LE.1d-6)*(-1d6)
cres = cres - (ctres.LE.1d-6)*(-1d6)

IF (sigma.EQ.1d0) THEN
    utilc = DLOG(MAX(cres,1d-8))
ELSE 
    utilc = MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma) 
ENDIF

IF (sigmag.EQ.1d0) THEN 
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag) 
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n +penalty 

objective_u_fe = wres

if (choice_glob.EQ.1) THEN
    gn_global1 = gnres
    x_global1  = hres
    ct_global1 = ctres
ENDIF

IF (try.EQ.1) THEN
    IF (objective_u_fe.GT.objective_u_feold) THEN
        objective_u_feold= objective_u_fe
        if (choice_glob.EQ.1) THEN
            gn_global1old = gn_global1    
            x_global1old  = x_global1 
            ct_global1old = ct_global1 
        ENDIF
    ENDIF

    !case 2: binding tax rate
    !pn*gn > 0
    
    temp1=-aux0/tau_value-a0
    IF (temp1.GT.0d0) THEN
        maxc = ((((1d0-omegac)/omegac)/temp1)**(1d0/(1d0+mu)))*ctval        !b2
        minc = 1d0-in*tau_value                                             !b6
    ELSE
        temp2 = temp1*tau_value/(in*tau_value-1d0)
        maxc  = (((1d0-omegac)/omegac/temp2)**(1d0/(1d0+mu)))*ctval         !b3
        maxc  = MIN(maxc,1d0-in*tau_value)                                  !b5
        maxc3=maxc
        IF (wbar.EQ.0d0) THEN             
            minc  = 1d-5                                                    !b4
        ELSE
            minc  = 1d0-in*taxrate-(aux0+taxrate*a0)*alpha/wbar             !b4
        ENDIF  

    ENDIF    
    IF (wbar.EQ.0d0) THEN
        maxc2  = 1d0             
    ELSE                
        !alpha*pn >= wbar
        maxc2 = ((1d0-omegac)/omegac*alpha/wbar)**(1d0/(1d0+mu))*ctval      !b1
    ENDIF

    maxc2 = MIN(maxc,maxc2)
    maxc  = maxc2

    ftype = 402
    IF (i_num_star.EQ.0) old_value = AFunction(ftype,a0,aux0,ct_global,minc,tau_value)
    old_valuemax = 10d+33 

    index_opt = 0
    index_optmax = 0
    difc = (maxc-minc)/(REAL(num)-1d0)

    150 DO i_num=1+i_num_star,num
        ftype=402
        ctemp(i_num) = minc+(REAL(i_num)-1d0)*difc     
        new_value = AFunction(ftype,a0,aux0,ct_global,ctemp(i_num),tau_value)

        sign_new = SIGN(1d0,new_value)
        sign_old = SIGN(1d0,old_value)
        sign_two = SIGN(1d0,sign_new*sign_old)

        if (DABS(new_value) <= DABS(old_valuemax)) then
            index_optmax = i_num
            old_valuemax = new_value        
        end if      
        
        i_num_star = i_num
        if (sign_two.EQ.-1d0) then
            index_opt = i_num
            old_value = new_value
            GOTO 222
        ENDIF
        old_value = new_value
    ENDDO       

    222 IF (index_opt.EQ.0) THEN
        bound1=ctemp(MAX(index_optmax-1,1))
        bound2=ctemp(MIN(index_optmax+1,num))
    ELSE
        bound1=ctemp(MAX(index_opt-1,1))
        bound2=ctemp(MIN(index_opt,num))        
    ENDIF

    IF (i_num_star.EQ.num)   try = 2
    
    cnval=BrentRoots(ftype,a0,aux0,ct_global,tau_value,bound1,bound2,1d-10,1000,fval,niter2,error)  
    gnval = 1d0-cnval
   
    IF (error.EQ.1) THEN 
        gnval = -10d0
        IF (i_num_star.LT.num)  go to 150
        GOTO 2046
    ELSEIF (error.EQ.2) THEN
        PRINT*,'Problem FE: More Iterations!'
        PAUSE
        STOP
    ENDIF

    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    wval   = alpha*p0
    tauval = tau_value*(a0+in*p0)
    
    GO TO 2046

ELSE
    IF (objective_u_feold.GT.objective_u_fe) THEN
        objective_u_fe= objective_u_feold
        if (choice_glob.EQ.1) THEN
            gn_global1 = gn_global1old    
            x_global1  = x_global1old 
            ct_global1 = ct_global1old 
        ENDIF
    ENDIF
ENDIF

end FUNCTION objective_u_fe
!***************************************************************************
REAL(8) FUNCTION objective_u_unemp(x,a_initial2,b_initial2,i_default_global2,tau_value)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::x,a_initial2,b_initial2,tau_value
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::pres,hres,gnres,cnres,ctres,cres,a0,b0,aux0
REAL(8)::wres,taures,ares,wwres,utilc,utilg,gnval,hval
REAL(8)::lambda1,pp1,bpval,wval,tauval
REAL(8)::p0,m0,h0,q,b,x_global2old,objective_u_unempold 
INTEGER,PARAMETER::num=100
REAL(8)::ctemp(num),aahat,aa1    
REAL(8)::ctval,cnval,cval,q_fun,penalty,fval,gn_global2old,ct_global,ct_global2old
REAL(8)::aux_lhs,balt,qalt,aux0alt,minc,maxc,temp1,taxbar,netbc,bound1,bound2,new_value,old_value,old_valuemax
INTEGER(4)::i,jj,iax,feas ,error,niter2,try,ftype,i_num,index_opt,i_num_star,i_num_star2,index_optmax,root
real(8)::sign_new,sign_old,sign_two,ctemp2,difc
LOGICAL::flag00
EXTERNAL q_fun

objective_u_unempold = -1d10
i_num_star2 = 0
root=0

b0 = b_initial2
a0 = DEXP(a_initial2)
b = x

q  = q_fun(b,a_initial2)

aux0 = -q*(b-(1d0-lambda)*b0) + payoff*b0            !p*gn - T
aux0 = aux0*(1d0-defaultgrid(i_default_global2))
aux_lhs   = aux0+tau_value*(a0+in*wbar/alpha)                          !always smaller than cT
!aux_lhs < 0 follows from  pn*gn+T = aux0+tax < 0 
!                          alpha*p*h**alpha=wbar*h < wbar

IF (choice_nowbar.EQ.1) aux_lhs = 100d0     !no aux_lhs condition for wbar=0

ctval  = a0+aux0
ct_global = ctval

IF (choice_fe.EQ.1) THEN
    hres = 0d0
    gnres= 0d0
    wres = huge_n
    flag00=(hres.LE.0d0)
    try  = 1
    GO TO 2033
ELSEIF (choice_nowbar.EQ.1) THEN
    objective_u_unemp= huge_n
    if (choice_glob.EQ.1) THEN
        gn_global2= 0d0    
        x_global2 = 0d0
        ct_global2= 0d0
    ENDIF
    RETURN
ELSEIF  ((aux_lhs.LE.0d0).AND.(tau_value.LT.1d0)) THEN
    balt = b+1d-6
    q    = q_fun(balt,a_initial2) 
    aux0alt = -q*(balt-(1d0-lambda)*b0) + payoff*b0            
    aux0alt = aux0alt*(1d0-defaultgrid(i_default_global2))
    objective_u_unemp= huge_n+1d3*(aux0alt-aux0)           
    if (choice_glob.EQ.1) THEN
        gn_global2 = 0d0    
        x_global2  = 0d0
        ct_global2 = 0d0
    ENDIF
    RETURN    
ENDIF

!case 1: h<1 and nonbinding tax cap
try = 1

ftype = 101
minc  = 1d-3
!Alternative: minc = (alpha/wbar*(1d0-omegac)/omegac)**(1d0/(1d0+mu))*ctval                       !b8
maxc  = (wbar/alpha*omegac/(1d0-omegac)/(ctval**(1d0+mu)))**(-alpha/(mu*(alpha-1d0)+1d0+mu))      !b7
!IF (choice_nowbar.EQ.1) maxc = 1d0

old_value = AFunction(ftype,a0,aux0,ct_global,minc,tau_value)
old_valuemax = 10d+33 

i_num_star=0
index_opt = 0
index_optmax = 0
difc = (maxc-minc)/(REAL(num)-1d0)

454 DO i_num=1+i_num_star,num
    ftype=101
    ctemp(i_num) = minc+(REAL(i_num)-1d0)*difc     
    new_value = AFunction(ftype,a0,aux0,ct_global,ctemp(i_num),tau_value)

    sign_new = SIGN(1d0,new_value)
    sign_old = SIGN(1d0,old_value)
    sign_two = SIGN(1d0,sign_new*sign_old)

    if (DABS(new_value) <= DABS(old_valuemax)) then
        index_optmax = i_num
        old_valuemax = new_value        
    end if      
    
    i_num_star = i_num
    if (sign_two.EQ.-1d0) then
        index_opt = i_num
        old_value = new_value
        GOTO 555
    ENDIF
    old_value = new_value
ENDDO       

555 IF (index_opt.EQ.0) THEN
    bound1=ctemp(MAX(index_optmax-1,1))
    bound2=ctemp(MIN(index_optmax+1,num))
ELSE
    bound1=ctemp(MAX(index_opt-1,1))
    bound2=ctemp(MIN(index_opt,num))        
ENDIF

cnval=BrentRoots(ftype,a0,aux0,ct_global,tau_value,bound1,bound2,1d-10,1000,fval,niter2,error)  

IF ((error.EQ.1)) THEN 
    cnval = -10d0
    IF (i_num_star.LT.num)  go to 454
ELSEIF (error.EQ.2) THEN
    PRINT*,'Problem: More Iterations!'
    PAUSE
    STOP
ENDIF

p0 = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
hval = (wbar/alpha/p0)**(1d0/(alpha-1d0))

1013 gnval  = hval**alpha-cnval
tauval = p0*gnval-aux0   

1015 taxbar = tau_value*(a0+in*p0*hval**alpha)
    
IF (chi.EQ.0d0) THEN
    flag00 = ((gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0)) 
ELSE
    flag00 = ((gnval.LT.0d0).OR.(gnval.GE.hval**alpha).OR.(p0.LE.0d0).OR.(ctval.LE.0d0).OR.(hval.GT.1d0)) 
ENDIF

IF (.NOT.flag00) THEN
    pres  = p0
    hres  = hval
    gnres = gnval
    ctres = ctval
    taures= tauval
    wwres = wbar
    feas  = 1d0
ELSE
    ctres = 0d0 
    feas  = 0d0
    hres  = 0d0
    gnres = 0d0
    taures= tauval
    pres  = p0
    wwres = wbar    
ENDIF
  
!negative penalty
penalty = 1d3*((pres.LE.0d0)*(0d0-pres)+(ctres.LE.0d0)*(0d0-ctres))
!without lumpsum taxes
IF (tau_value.LE.1d0) penalty = penalty+1d3*((taures.GT.taxbar)*(taures-taxbar))
penalty = penalty+1d9*((b.LT.bmin)*(bmin-b))+1d2*((b.LT.bmin))
penalty = penalty +1d3*((gnres.GE.hres**alpha)*(gnres-hres**alpha)+(gnres.LT.0d0)*(0d0-gnres))
penalty = penalty +1d3*((hres.GT.1d0)*(hres-1d0)+(hres.LE.0d0)*(0d0-hres))

cnres = MAX(1d-8,hres**alpha-gnres)

IF (mu .NE. 0d0) THEN
    cres = (omegac*ctres**(-mu)+ (1d0-omegac)*cnres**(-mu))**(-1d0/mu)
ELSE
    cres = (ctres**(omegac))*(cnres**(1d0-omegac))
ENDIF

IF ((cnres.LE.1d-6).OR.(ctres.LE.1d-6)) cres = 0d0

IF (alphac.EQ.1d0) THEN
    aahat  = 1d0
    aa1    = 1d0
ELSE
    aahat  = (hres+(1d0-hres)*alphac)**(-1d0)
    aa1    = hres+(1d0-hres)*(alphac**(1d0-sigma))
ENDIF

    
IF (sigma.EQ.1d0) THEN
    utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(cres,1d-8))
ELSE    
    utilc = (aahat**(1d0-sigma))*aa1*MAX(cres,1d-8)**(1d0-sigma)/(1d0-sigma)
ENDIF

IF (sigmag.EQ.1d0) THEN
    utilg = DLOG(MAX(gnres,1d-8))
ELSE
    utilg = MAX(gnres,1d-8)**(1d0-sigmag)/(1d0-sigmag)
ENDIF

wres = (1d0-chi)*utilc+chi*utilg 
wres = feas*wres + (1d0-feas)*huge_n + penalty 

objective_u_unemp = wres
if (choice_glob.EQ.1) THEN
    gn_global2 = gnres
    x_global2 = hres
    ct_global2 = ctres
ENDIF
        
IF (try.EQ.1) THEN
    IF (root.EQ.0) THEN
        x_global2old  = x_global2
        gn_global2old = gn_global2
        ct_global2old = ct_global2
    
        objective_u_unempold = objective_u_unemp
        root=1
    ELSEIF (root.GT.0) THEN
        IF (objective_u_unemp.GT.objective_u_unempold) THEN
            objective_u_unempold= objective_u_unemp
        
            if (choice_glob.EQ.1) THEN
                gn_global2old = gn_global2    
                x_global2old  = x_global2 
                ct_global2old = ct_global2
            ENDIF
        ENDIF
        root=root+1
    ENDIF
    
    IF (i_num_star.LT.num)     GO TO 454
ENDIF

2033 IF (try.EQ.1) THEN
    
    !case2: h=1 and nonbinding cap         
    try = 2
    IF (choice_fe.EQ.1) try = 3

    hval  = 1d0
    p0    = wbar/alpha
    ctval = ct_global
    cnval = (((1d0-omegac)/omegac/p0)**(1d0/(1d0+mu)))*ctval

    IF (tau_value.GT.1d0) THEN
        IF (chi.EQ.0d0) THEN
            try=4
        ELSE
            try=5
        ENDIF
    ENDIF

    GO TO 1013
    
ELSEIF ((try.GT.1).AND.(try.LT.4)) THEN  
      
    IF (objective_u_unemp.GT.objective_u_unempold) THEN
        objective_u_unempold= objective_u_unemp
        
        if (choice_glob.EQ.1) THEN
            gn_global2old = gn_global2    
            x_global2old  = x_global2 
            ct_global2old = ct_global2
        ENDIF
    ENDIF

    IF (try.EQ.2) THEN

        !case 3: h<1 and binding tax cap
        IF (tau_value.GT.1d0) THEN
            try =100
        ELSE
            IF (i_num_star2.EQ.num) try = 3
        ENDIF
        
        maxc = ((1d0-omegac)/omegac*alpha/wbar*(ctval**(1d0+mu)))**(alpha/(1d0+alpha*mu))       !b7        
        !maxc=MIN(1d0-1d-6,maxc)
        
        !h < 1
        minc = (alpha/wbar*(1d0-omegac)/omegac)**(1d0/(1d0+mu))*ctval                           !b8
        !minc = MAX(1d-5,minc)
 !       IF (choice_nowbar.EQ.1) THEN
 !           minc = 1d-3
 !           maxc = 1d0
 !        ENDIF
        
        ftype = 302
        IF (i_num_star2.EQ.0)  old_value = AFunction(ftype,a0,aux0,ct_global,minc,tau_value)
        old_valuemax = 10d+33 

        i_num_star2=i_num_star2
        index_opt = 0
        index_optmax = 0
        difc = (maxc-minc)/(REAL(num)-1d0)

        343 DO i_num=1+i_num_star2,num
            ftype = 302
            ctemp(i_num) = minc+(REAL(i_num)-1d0)*difc     
            new_value = AFunction(ftype,a0,aux0,ct_global,ctemp(i_num),tau_value)

            sign_new = SIGN(1d0,new_value)
            sign_old = SIGN(1d0,old_value)
            sign_two = SIGN(1d0,sign_new*sign_old)

            if (DABS(new_value) <= DABS(old_valuemax)) then
                index_optmax = i_num
                old_valuemax = new_value        
            end if      
        
            i_num_star2 = i_num
            if (sign_two.EQ.-1d0) then
                index_opt = i_num
                old_value = new_value
            GOTO 444
            ENDIF
            old_value = new_value
        ENDDO       
 
        444 IF (index_opt.EQ.0) THEN
            bound1=ctemp(MAX(index_optmax-1,1))
            bound2=ctemp(MIN(index_optmax+1,num))
        ELSE
            bound1=ctemp(MAX(index_opt-1,1))
            bound2=ctemp(MIN(index_opt,num))        
        ENDIF

        cnval=BrentRoots(ftype,a0,aux0,ct_global,tau_value,bound1,bound2,1d-10,1000,fval,niter2,error)  

        IF ((error.EQ.1)) THEN 
            cnval = -10d0
            gnval = -10d0

            IF (i_num_star2.LT.num)  go to 343
            GOTO 1015
        ELSEIF (error.EQ.2) THEN
            PRINT*,'Problem: More Iterations!'
            PAUSE
            STOP
        ENDIF
        
        p0    = (1d0-omegac)/omegac*(ct_global/cnval)**(1d0+mu)
        hval  = (alpha/wbar*p0)**(1d0/(1d0-alpha))

        gnval = hval**alpha-cnval
        tauval= tau_value*(a0+in*p0*hval**alpha)
        
        GO TO 1015     
    ELSEIF (try.EQ.3) THEN         
        !case 4: h=1 and binding tax cap
        IF (chi.EQ.0d0) THEN
            try = 4
        ELSE
            try=5
        ENDIF
                
        hval  = 1d0
        p0    = wbar/alpha
        temp1 = (p0*omegac/(1d0-omegac))**(1d0/(1d0+mu))

        cnval = (a0-tau_value*(a0+in*p0)+p0)/(temp1+p0)
        ctval = temp1*cnval
        gnval = 1d0-cnval
                tauval=p0*gnval-aux0
        taxbar= taxrate*(a0+in*p0)
        
        netbc = tauval-taxbar 

        IF (DABS(netbc).GT.1d-6) tauval=tauval+10d0     !need to impose govt bc with equality

        GO TO 1015
    ENDIF
ELSEIF (try.EQ.4) THEN
    IF (objective_u_unemp.GT.objective_u_unempold) THEN
        objective_u_unempold= objective_u_unemp
        
        if (choice_glob.EQ.1) THEN
            gn_global2old = gn_global2    
            x_global2old  = x_global2 
            ct_global2old = ct_global2
        ENDIF
    ENDIF
    ctval=ct_global
    hval = (alpha/wbar*(1d0-omegac)/omegac*(ctval**(1d0+mu)))**(1d0/(1d0+mu*alpha))

    IF ((hval.LE.0d0).OR.(hval.GT.1d0)) THEN
        hval=0d0
        cnval=0d0
        p0=0d0
    ELSE
        cnval=hval**alpha
        p0 = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    ENDIF

    try=5
    gnval  = 0d0
    tauval = -aux0 

    GO TO 1015
ELSE
    IF (objective_u_unempold.GT.objective_u_unemp) THEN
        objective_u_unemp= objective_u_unempold

        if (choice_glob.EQ.1) THEN
            gn_global2 = gn_global2old    
            x_global2  = x_global2old 
            ct_global2 = ct_global2old 
        ENDIF
    ENDIF
ENDIF


end function objective_u_unemp
!******************************************************************************************
REAL(8) FUNCTION objective_u_tax(b,a_initial2,b_initial2,i_default_global2,tau_value2)

USE GLOBAL

IMPLICIT NONE

REAL(8),INTENT(IN)::b,a_initial2,b_initial2,tau_value2
INTEGER,INTENT(IN)::i_default_global2
REAL(8)::value1,value2,bb,objective_u_fe,objective_u_unemp
EXTERNAL::objective_u_fe,objective_u_unemp

bb = b*(1d0-defaultgrid(i_default_global2))

value1= objective_u_fe(bb,a_initial2,b_initial2,i_default_global2,tau_value2)
value2= objective_u_unemp(bb,a_initial2,b_initial2,i_default_global2,tau_value2)

objective_u_tax = MAX(value1,value2)

if (choice_glob.EQ.1) THEN
    IF (value1.GE.value2) THEN
        reg_global = 1
        x_global    = x_global1
        gn_global    = gn_global1
        ctrad_global = ct_global1  
    ELSE
        reg_global = 2
        x_global     = x_global2
        gn_global    = gn_global2
        ctrad_global = ct_global2  
    ENDIF    
ENDIF

END FUNCTION objective_u_tax
!**************************************************************************************
!Compute the objective function in the Belman equation
DOUBLE PRECISION function objective_fun(b,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

IMPLICIT NONE

DOUBLE PRECISION,INTENT(IN)::b,a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
INTEGER :: other_type, i,j,h, t, d, NINTV, index_opt
REAL(8)::u_fun,q_fun,acum,exp_v_next,value_next,a_next,a_fun,scalar,DCSVAL,bb
REAL(8)::a_threshold,a_max1,a_min1,a_next1,acum2,acum1,g,q,DNORDF,v0_fun,vaut_fun,value_next_tirar
REAL(8)::acum_int,acum_int_lo,acum_int_hi,borrowing, interpolate,cdf,cdf1,cdf_threshold, DNORIN,objective_u_tax
REAL(8)::exp_v_next_aut,acum3
EXTERNAL u_fun,q_fun,DNORDF,a_fun,interpolate,v0_fun,vaut_fun,DCSVAL,DNORIN,udef,objective_u_tax

REAL(8)::value1,value2 ,u_value,ucost,udef,old_value,new_value, tau_value, tau_value2

IF (psi1.LT.500d0) THEN
    a_threshold   = a_fun(b,a_initial2)
ELSE
    a_threshold   = -1d4
ENDIF

cdf_threshold = DNORDF((a_threshold - rhoa*a_initial2 - (1-rhoa)* mua)/sigmaa)
bb = b*(1d0-defaultgrid(i_default_global2))

exp_v_next = 0d0
exp_v_next_aut = 0d0
scalar = 1d0 / SQRT(pi_number)

a_min1 = (1-rhoa)*mua + rhoa* a_initial2 - 4*sigmaa
a_max1 = (1-rhoa)*mua + rhoa* a_initial2 + 4*sigmaa

NINTV = ncdf - 1
acum1=0d+0
acum2=0d+0

if (a_threshold < a_min1) then

    do i=1,nquad_v
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = v0_fun(bb,a_next)
        acum1 = acum1 + 0.5*(cdfmax-cdfmin)*quad_w_v(i) * value_next
    end do

elseif(a_threshold > a_max1) then

    do i=1,neps
        cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
        a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
        value_next = vaut_fun(a_next)
        acum1 = acum1 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
    end do

else
    do i=1,neps

    cdf = 0.5*(quad_x_v(i) + 1)*(cdf_threshold - cdfmin) + cdfmin
    a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
    value_next = vaut_fun(a_next)
    acum1 = acum1 + 0.5*(cdf_threshold - cdfmin)*quad_w_v(i) * value_next

!        compute integral for a > a_threshold
    cdf1 = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdf_threshold) + cdf_threshold
    a_next1 = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf1) * sigmaa
    value_next = v0_fun(bb,a_next1)
    acum2 = acum2 + 0.5*(cdfmax - cdf_threshold)*quad_w_v(i) * value_next

    end do
end if

exp_v_next = (acum1+acum2)/(cdfmax - cdfmin)

IF (i_default_global2.EQ.2) THEN
    acum3=0d+0
        do i=1,neps
            cdf = 0.5*(quad_x_v(i) + 1)*(cdfmax - cdfmin) + cdfmin
            a_next  = (1-rhoa)*mua + rhoa* a_initial2 + DNORIN(cdf) * sigmaa
            value_next = vaut_fun(a_next)
            acum3 = acum3 + 0.5*(cdfmax - cdfmin)*quad_w_v(i) * value_next
        end do
    
        exp_v_next_aut = acum3/(cdfmax - cdfmin)
        exp_v_next = phi*exp_v_next + (1d0-phi)*exp_v_next_aut
ENDIF

ucost = defaultgrid(i_default_global2)*udef(a_initial2)

old_value=-10d33

DO i=1,ntau
    tau_value2 = taugrid(i)

   new_value = objective_u_tax(b,a_initial2,b_initial2,i_default_global2,tau_value2) - taucost*(tau_value2-taxrate)**2d0
   if (new_value >= old_value) then
      index_opt = i
      old_value = new_value
  end if
end do

tau_value = taugrid(index_opt)
tau_global = tau_value

u_value = objective_u_tax(b,a_initial2,b_initial2,i_default_global2,tau_value) - taucost*(tau_value-taxrate)**2d0

objective_fun = -(u_value - ucost) - beta * exp_v_next

end FUNCTION objective_fun
!********************************************************************************************
subroutine optimize(b_next,v_value,a_initial2,b_initial2,i_default_global2)

USE GLOBAL

implicit none

REAL(8),INTENT(IN)::a_initial2,b_initial2
INTEGER,INTENT(IN)::i_default_global2
integer :: MAXFN, i,t, d, index_opt, index_vector(1), num, num_tirar
parameter (num=100, num_tirar = 500)      
DOUBLE PRECISION :: b_next, v_value, objective_fun, new_value, q_fun, old_value, b, u_fun,q,vector(nb),&
                    b_next_grid(num), STEP, BOUND, XACC, b_next_guess
real(8)::b_next_inf,b_next_sup 
EXTERNAL objective_fun, q_fun, u_fun, DUVMIF

REAL(8)::objective_u_unemp,objective_u_fe,value1,value2
EXTERNAL::objective_u_unemp,objective_u_fe

IF (psi1.GT.500d0) THEN
    b_next_inf = MAX(bmin,b_initial2-0.2d0) 
    b_next_sup = MIN(bmax,b_initial2+0.2d0)
ELSE
    b_next_inf = bmin 
    b_next_sup = bmax
ENDIF

old_value = 10d+33
do i=1,num 

   b_next_grid(i) = b_next_inf + (b_next_sup - b_next_inf) * (i-1)/ (num - 1)
   new_value = objective_fun2(b_next_grid(i))
   
   if (new_value <= old_value) then
      index_opt = i
      b_next_guess = b_next_grid(i)
      old_value = new_value
  end if
end do

STEP = (b_next_grid(2) - b_next_grid(1))*0.5
BOUND = 10*EXP(amax)
XACC = 1d-6
MAXFN = 1000
   
IF (old_value.GE.1d5) THEN
    b_next = 0d0
    v_value = -old_value       
    else
        call DUVMIF(objective_fun2, b_next_guess, STEP, BOUND, XACC, MAXFN, b_next)
        v_value = -objective_fun2(b_next)
        !value1 = objective_u_fe(b_next,a_initial2,b_initial2,i_default_global2)
        !value2 = objective_u_unemp(b_next,a_initial2,b_initial2,i_default_global2)     
end if

CONTAINS

REAL(8) FUNCTION objective_fun2(x)

REAL(8),INTENT(IN)::x

objective_fun2 = objective_fun(x,a_initial2,b_initial2,i_default_global2)

END FUNCTION
100 end subroutine
!**********************************************************************************
subroutine iterate

USE GLOBAL
USE OMP_LIB

IMPLICIT NONE

INTEGER :: d

REAL(8)::a_initial,b_initial
INTEGER::i_default_global
REAL(8)::a_valor,b_valor,b_valor2,v_valor,convergence,criteria,deviation,q_fun,b_next,bmin_old
REAL(8)::b0_next,b_next1(na),b,w_value,dev_q_paid,dev_q_scheme,dev_v
REAL(8)::v0_value,v1_value,q,interpolate,b_next0,b_grid_local(nb),bmin_local
REAL(8)::FDATA(nb),DLEFT,DRIGHT,BREAK_GRID(nb),CSCOEF(4,nb),DCSVAL,epsilongrid(ncdf), dev_q_nodef

DOUBLE PRECISION, DIMENSION(nb,na)::v0_matrix_new, vaut_matrix_new, q_matrix_new, default_decision_new,q_scheme, q_scheme_new
DOUBLE PRECISION, DIMENSION(nb,na)::v_matrix_new,dev_matrix,q_matrix_nodef_new
REAL(8)::yTmax,value1
INTEGER i_b, i_a, i, j, i_b_zero, indices(1:2), indice_b, counter, ILEFT, IRIGHT,NINTV,DNORIN,new_index_b,index_b !i_def_opt
REAL(8)::v0_fun,x0,x1(na),b_fun,bmin_feas,objective_u_fe,objective_u_unemp,objective_u_tax,objective_fun
REAL(8)::min_m,natlevel,value2,gntemp(na),cttemp(na),tautemp(na),ct_matrix_new(nb,na),tau_matrix_new(nb,na)
EXTERNAL::objective_u_fe,objective_u_unemp,objective_u_tax,objective_fun
REAL(8)::new_value
INTEGER::index_bfeas(na),index_opt,ii_b,nbaux
CHARACTER(len=29)::filename

criteria = tol  
convergence = -1
counter = 0
bmin_old = bmin
vaut_matrix_new = vaut_matrix
v0_matrix_new = v0_matrix
v_matrix_new = v_matrix
q_scheme = 0d+0

do WHILE((convergence<0).AND.(counter.LE.maxiter))
    index_bfeas4 = 0d0
    
    !to compute spline for feasible debt gridpoints
    do i_a = 1,na
        new_value=-1d5
        i_b=0
        index_opt = 0
        
        do while (new_value <= -1d3)
           index_opt = i_b
           IF (i_b.EQ.nb) EXIT
           i_b=i_b+1
           new_value = v0_matrix(i_b,i_a)
        end do
        index_bfeas4(i_a)=index_opt     
    ENDDO

    index_bfeas2=MAXVAL(index_bfeas4)

      deviation = 0d+0
      counter = counter + 1
      dev_v   = 0d+0
      dev_q_scheme = 0d+0
      dev_q_paid   = 0d+0
      dev_q_nodef  = 0d+0    
         
     do i_a = 1,na         
         IF (index_bfeas4(i_a)+1.LE.nb-1) THEN      !need at least to points for splines
             ii_b=index_bfeas4(i_a)+1
             nbaux =nb-ii_b+1
             FDATA(1:nbaux) = v0_matrix(ii_b:nb,i_a)
	         ILEFT  = 0
	         IRIGHT = 0
	         DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(ii_b+1) - bgrid(ii_b))
	         DRIGHT = (FDATA(nbaux)-FDATA(nbaux-1)) / (bgrid(nbaux) - bgrid(nbaux-1))
	         call  DCSDEC (nbaux, bgrid(ii_b:nb), FDATA(1:nbaux), ILEFT, DLEFT, IRIGHT, DRIGHT, break_matrix(i_a, ii_b:nb), coeff_matrix(i_a, :,ii_b:nb))
         ENDIF         
     ENDDO
         
        do i_a = 1,na
            FDATA = q_matrix_nodef(:,i_a)
	        ILEFT  = 0
	        IRIGHT = 0
            DLEFT = (FDATA(2)-FDATA(1)) / (bgrid(2) - bgrid(1))
	        DRIGHT = (FDATA(nb)-FDATA(nb-1)) / (bgrid(nb) - bgrid(nb-1))
	        call  DCSDEC (nb, bgrid, FDATA, ILEFT, DLEFT, IRIGHT, DRIGHT, BREAK_GRID, CSCOEF)
            break_matrix_q(i_a, :) = BREAK_GRID
            coeff_matrix_q(i_a, :,:) = CSCOEF
        end do

    IF (psi1.LE.500d0) THEN
    !! Compute the value of repayment
    !!!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(a_initial,b_initial,i_default_global,b_valor,v_valor)
    !!!$OMP DO COLLAPSE(1)
    LOOP_A1: do i_a = 1,na
                  a_initial = agrid(i_a)
                  b_initial = bgrid(1)
    
                  i_default_global=2 !Country defaults
                  b_valor = 0d0
                  v_valor= -objective_fun(b_valor,a_initial,b_initial,i_default_global)
                  
                  x1(i_a)     = x_global
                  gntemp(i_a) = gn_global
                  tautemp(i_a) = tau_global
                  cttemp(i_a) = ctrad_global
                  b_next1(i_a) = b_valor
                  vaut_matrix_new(:, i_a) = v_valor
            end do LOOP_A1
            !!!$OMP END DO
            !!!$OMP END PARALLEL
            ELSE
                  x1      = 0d0
                  gntemp  = 0d0
                  tautemp = 0d0
                  cttemp  = 0d0
                  b_next1 = 0d0
                  vaut_matrix_new = -1d10
            ENDIF
    
    !!!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(a_initial,b_initial,i_default_global,v_valor,b_next0)
    !!!$OMP DO COLLAPSE(2)
    LOOP_A2: do i_a = 1,na
    LOOP_B:  do i_b = 1,nb 
                  a_initial = agrid(i_a)
                  b_initial = bgrid(i_b)

                  i_default_global=1 !Country does not default
                  
                  IF (welfare_case.EQ.1) THEN
                        b_next0 = (1d0-lambda)*b_initial
                        v_valor= -objective_fun(b_next0,a_initial,b_initial,i_default_global)
                  ELSE
                    call optimize(b_next0, v_valor,a_initial,b_initial,i_default_global)
                  ENDIF
                  v0_matrix_new(i_b, i_a) = v_valor

                  q_matrix_nodef_new(i_b, i_a) = q_fun(b_next0, a_initial)
                  x0 = x_global
                  ct_matrix_new(i_b, i_a) = ctrad_global    
                  
                  if (vaut_matrix_new(i_b, i_a) > v0_matrix_new(i_b, i_a)) THEN 
                     b_next_matrix(i_b, i_a) = 0d0 
                     x_matrix_new(i_b, i_a)  = x1(i_a)                     
                     gn_matrix_new(i_b, i_a) = gntemp(i_a)     
                     tau_matrix_new(i_b, i_a)= gntemp(i_a)     
                     default_decision_new(i_b,i_a) = 2
                     v_matrix_new(i_b, i_a) = vaut_matrix_new(i_b, i_a)
                  else
                        b_next_matrix(i_b, i_a) = b_next0                     
                        gn_matrix_new(i_b, i_a) = gn_global   
                        x_matrix_new(i_b, i_a)  = x0                     
                        tau_matrix_new(i_b, i_a)= tau_global    
                        
                        default_decision_new(i_b, i_a) = 1
                        v_matrix_new(i_b, i_a) = v0_matrix_new(i_b, i_a)
                  END if
                  q_scheme_new(i_b, i_a) =  q_fun(b_initial,a_initial)
                  
                   q_matrix_new(i_b, i_a) = q_fun(b_next_matrix(i_b,i_a),a_initial)               

                   dev_v =        max(ABS(v_matrix_new(i_b, i_a) - v_matrix(i_b, i_a)), dev_v)
                   dev_q_paid =   MAX(ABS(q_matrix_new(i_b, i_a) - q_matrix(i_b, i_a)), dev_q_paid)
                   dev_q_scheme = MAX(ABS(q_scheme_new(i_b, i_a) - q_scheme(i_b, i_a)), dev_q_scheme)
                   dev_q_nodef =   MAX(ABS(q_matrix_nodef_new(i_b, i_a) - q_matrix_nodef(i_b, i_a)), dev_q_nodef)

                   deviation = MAX(deviation, MAX(dev_v, dev_q_scheme))
                   dev_matrix(i_b, i_a) = (q_matrix_new( i_b, i_a) - q_matrix( i_b, i_a))
                end do LOOP_B
     end do LOOP_A2
    !!!$OMP END DO
    !!!$OMP END PARALLEL
                
       bmin_local = bmin


 WRITE(*, '(I5, X, F15.8, X, F12.8, X, F12.8, X, F12.8)') counter, dev_v,dev_q_scheme,MINVAL(-b_next_matrix*(v0_matrix.GE.-500d0)),MAXVAL(b_next_matrix)

bmin_local = bmin_local 
    
!save results of the current iteration
IF (save_results.EQ.1) THEN
      open (10, FILE='graphs/v.txt',STATUS='replace')
      open (11, FILE='graphs/default.txt',STATUS='replace')
      open (12, FILE='graphs/q.txt', STATUS='replace')
      open (13, FILE='graphs/b_next.txt',STATUS='replace')
      open (14, FILE='graphs/devq.txt',STATUS='replace')
      open (15, FILE='graphs/dev_iter.txt',ACCESS='append')
      open (16, FILE='graphs/q_paid.txt', STATUS='replace')
      open (110, FILE='graphs/b_grid.txt',STATUS='replace')
      open (113, FILE='graphs/bounds.txt',STATUS='replace')
      open (114, FILE='graphs/x.txt',STATUS='replace')
      open (115, FILE='graphs/gn.txt',STATUS='replace')
      open (116, FILE='graphs/ctt.txt',STATUS='replace')
      open (117, FILE='graphs/tau.txt',STATUS='replace')
      
      WRITE(113, '(F15.11, X, F15.11)') bmin, bmax
      WRITE(113, '(F15.11, X, F15.11)') amin, amax

      WRITE(15,'(I3,X,F15.10, X, F15.10)') counter,dev_v ,dev_q_paid

         do i_b = 1,nb
            WRITE(110, '(F15.11)') bgrid(i_b)
            do i_a = 1,na
                     a_initial = agrid(i_a)
                     b_initial = bgrid(i_b)

                     WRITE(10, '(F18.10, X, F18.10, X, F18.10)') MAX(-1000d0,v_matrix_new(i_b, i_a)), &
                     MAX(-1000d0,v0_matrix_new(i_b, i_a)), MAX(-1000d0,vaut_matrix_new(i_b, i_a))
                     WRITE(11, '(F6.2)') defaultgrid(default_decision(i_b, i_a))

                     WRITE(12, '(F15.11)') q_fun(bgrid(i_b), agrid(i_a))
                     WRITE(13, '(F15.11)') b_next_matrix(i_b, i_a)
                     WRITE(14, '(F15.11)') dev_matrix(i_b, i_a)
                     WRITE(16, '(F15.11, X, F15.11)') q_matrix_new(i_b, i_a), q_matrix_nodef_new(i_b, i_a)
                     WRITE(114, '(F15.11)') x_matrix_new(i_b, i_a)
                     WRITE(115, '(F15.11)') gn_matrix_new(i_b, i_a)                     
                     WRITE(116, '(F15.11)') ct_matrix_new(i_b, i_a)                     
                     WRITE(117, '(F15.11)') tau_matrix_new(i_b, i_a)                     
            end do
         end do

 CLOSE(10)
 CLOSE(11)
 CLOSE(12)
 CLOSE(13)
 CLOSE(14)
 CLOSE(15)
 CLOSE(16)
 CLOSE(110)
 CLOSE(113)
 CLOSE(114)
 CLOSE(115)
 CLOSE(116)
 CLOSE(117)

ENDIF



!UPDATE VALUES OF MATRICES
   default_decision = default_decision_new
   if (deviation < criteria) then
      convergence =1
   end if
    
   IF (choice_bound.EQ.0) bmin_local=bmin

   if (bmin_local == bmin) then
       q_matrix = q_matrix_new
       q_scheme = q_scheme_new
       q_matrix_nodef = q_matrix_nodef_new
       vaut_matrix = weightv*vaut_matrix_new+(1-weightv)*vaut_matrix
       v0_matrix = weightv*v0_matrix_new+(1d0-weightv)*v0_matrix
       v_matrix = max(v0_matrix, vaut_matrix)
    else

   !UPDATE LOCAL GRID FOR FOR b
   do i=1,nb
       b_grid_local(i) = bmin_local + (bmax - bmin_local) * (i-1) / (nb - 1)
   end do

   do j=1,na
       vaut_matrix(:,j) = interpolate(b_grid_local(1), agrid(j), vaut_matrix_new)
       do i=1,nb
           v0_matrix(i,j) = interpolate(b_grid_local(i), agrid(j), v0_matrix_new)
           v_matrix(i,j) = max(v0_matrix(i,j), vaut_matrix(i,j))
       end do
   end do
   do i_b=1,nb
       do i_a =1,na
           q_matrix(i_b, i_a) = interpolate(b_grid_local(i_b), agrid(i_a), q_matrix_new)
           q_scheme(i_b, i_a) = interpolate(b_grid_local(i_b), agrid(i_a), q_scheme_new)
           q_matrix_nodef(i_b, i_a) = interpolate(b_grid_local(i_b), agrid(i_a), q_matrix_nodef_new)
       end do
   end do

   bmin_old = bmin
   bmin = bmin_local
   call compute_bgrid
    end if

    IF (counter.EQ.maxiter) EXIT

end do

counter_global = counter
 
PRINT*, 'Model solved'

end subroutine iterate
!************************************************************************************************************
SUBROUTINE pol_fcns
USE GLOBAL

IMPLICIT NONE

REAL(8)::a_initial,b_initial,a0,aux0
INTEGER::i_default_global
REAL(8)::a,b
INTEGER::i_a,i_b,i_bp,i,j,i0(1)
REAL(8)::wmax,p0max,h0max,gn0max,ct0max,tau0max,ww0max,q_fun,objective_fun,interpolate,v_valor,b_next
REAL(8)::ctval,cnval,tauval,p0,q,ca1(na2,nb2),q1_nodef(na2,nb2)
EXTERNAL::q_fun,objective_fun,interpolate

do i=1,nb2
   bgrid2(i) = bmin + (bmax - bmin) * REAL((i-1d0) / (nb2 - 1d0))
end do

do i=1,na2
   agrid2(i) = amin + (amax - amin) * REAL((i-1d0) / (na2 - 1d0))
end do

i0 = LOCATE2(bgrid2,0d0)
i_bzero2 = i0(1)
bgrid2(i_bzero2)=0d0
choice_glob=1

OPEN(1,FILE='output//bgrid2.txt')           
DO i = 1,nb2
    WRITE(1,*) bgrid2(i)
ENDDO
CLOSE(1)

OPEN(1,FILE='output//agrid2.txt')           
DO i = 1,na2
    WRITE(1,*) agrid2(i)
ENDDO    
CLOSE(1)

CALL BUILD_TRANSITION(na2,agrid2,rhoa,sigmaa,mua,prob)
   
DO i_a=1,na2

    a_initial = agrid2(i_a)
    a0 = DEXP(a_initial)

    b_initial = 0d0 
    i_default_global = 2  !default occurs
    b = 0d0

    v_valor = -objective_fun(b,a_initial,b_initial,i_default_global)
    vaut1(i_a)  = v_valor
    ctval  = ctrad_global 
    cnval  = x_global**alpha-gn_global
    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)
    tauval = tau_global
    
    paut1(i_a)  = p0
    haut1(i_a)  = x_global
    gnaut1(i_a) = gn_global
    ctaut1(i_a) = ctval
    
    tauaut1(i_a)= tauval*(a0+in*p0*x_global**alpha)
    tauvalaut1(i_a)=tauval
    waut1(i_a)  = alpha*p0*x_global**(alpha-1d0)
    
    DO i_b=1,nb2    
            
    i_default_global = 1
    b_initial = bgrid2(i_b)
    
    b = interpolate(b_initial, a_initial, b_next_matrix)
    i_bp = LOCATE2(bgrid2,b)
    
    IF (welfare_case.EQ.1) THEN
        b_next = (1d0-lambda)*b_initial
        v_valor= -objective_fun(b_next,a_initial,b_initial,i_default_global)
    ELSE
        call optimize(b_next, v_valor,a_initial,b_initial,i_default_global)
    ENDIF
    
    vcon1(i_a,i_b)  = v_valor
    bnext1(i_a,i_b) = b_next

    q = q_fun(b_next,a_initial)
    aux0 = -q*(b_next-(1d0-lambda)*b_initial ) + payoff*b_initial            !p*gn - T

    ctval  = ctrad_global
    cnval  = x_global**alpha-gn_global
    p0     = (1d0-omegac)/omegac*(ctval/cnval)**(1d0+mu)    
    bp1(i_a,i_b)    = LOCATE2(bgrid2,b_next)
    tauval = tau_global    
    q1_nodef(i_a,i_b) = q_fun(b_next,a_initial)
    p1(i_a,i_b)     = p0
    h1(i_a,i_b)     = x_global
    gn1(i_a,i_b)    = gn_global
    ct1(i_a,i_b)    = ctval
    q1(i_a,i_b)     = q_fun(b_initial,a_initial)
    ca1(i_a,i_b)    = aux0
    tauval1(i_a,i_b)= tauval
    tau1(i_a,i_b)   = tauval*(a0+in*p0*x_global**alpha)
    w1(i_a,i_b)     = alpha*p0*x_global**(alpha-1d0)                    

    def1(i_a,i_b) = 1d0
    IF (vaut1(i_a).GT.vcon1(i_a,i_b)) def1(i_a,i_b) = 0d0
ENDDO
ENDDO


OPEN(1,FILE='output//gvcon.txt')
OPEN(2,FILE='output//gp.txt')
OPEN(3,FILE='output//gh.txt')
OPEN(4,FILE='output//gibp.txt')
OPEN(5,FILE='output//gct.txt')
OPEN(6,FILE='output//gbnext.txt')
OPEN(7,FILE='output//ggn.txt')
OPEN(8,FILE='output//gtau.txt')    
OPEN(88,FILE='output//gtauval.txt')    
OPEN(9,FILE='output//gw.txt')    
OPEN(10,FILE='output//gq.txt')    
OPEN(11,FILE='output//gdef1.txt')    
OPEN(12,FILE='output//gca.txt')    
OPEN(13,FILE='output//gq_nodef.txt')    

DO i = 1,na2
DO j = 1,nb2
    WRITE(1,*) vcon1(i,j)
    WRITE(2,*) p1(i,j)
    WRITE(3,*) h1(i,j)
    WRITE(4,'(I7)') bp1(i,j)
    WRITE(6,*) bnext1(i,j)
    WRITE(5,*) ct1(i,j)
    WRITE(7,*) gn1(i,j)
    WRITE(8,*) tau1(i,j)
    WRITE(88,*) tauval1(i,j)
    WRITE(9,*) w1(i,j)
    WRITE(10,*) q1(i,j)
    WRITE(11,*) def1(i,j)
    WRITE(12,*) ca1(i,j)
    WRITE(13,*) q1_nodef(i,j)
ENDDO
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(4)
CLOSE(5)
CLOSE(6)
CLOSE(7)
CLOSE(8)
CLOSE(88)
CLOSE(9)
CLOSE(10)
CLOSE(11)
CLOSE(12)
CLOSE(13)

OPEN(1,FILE='output//gvaut.txt')
OPEN(2,FILE='output//gpaut.txt')
OPEN(3,FILE='output//ghaut.txt')
OPEN(5,FILE='output//gctaut.txt')
OPEN(7,FILE='output//ggnaut.txt')
OPEN(8,FILE='output//gtauaut.txt')    
OPEN(88,FILE='output//gtauvalaut.txt')    
OPEN(9,FILE='output//gwaut.txt')    
            
DO i = 1,na2
    WRITE(1,*) vaut1(i)
    WRITE(2,*) paut1(i)
    WRITE(3,*) haut1(i)
    WRITE(5,*) ctaut1(i)
    WRITE(7,*) gnaut1(i)
    WRITE(8,*) tauaut1(i)
    WRITE(88,*) tauvalaut1(i)
    WRITE(9,*) waut1(i)
ENDDO
    
CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(5)
CLOSE(7)
CLOSE(8)
CLOSE(88)
CLOSE(9)
    
END SUBROUTINE pol_fcns
!****************************************************************************************
SUBROUTINE MC_SIM(xvector)
! compute conditional moments

USE GLOBAL

IMPLICIT NONE

REAL(8),DIMENSION(6),INTENT(OUT)::xvector
REAL(8),DIMENSION(nt)::at,pt,rert,ht,ctt,gnt,bt,acct,deft,rt,wt,taut,tauvalt,bort,fmt,dtauvalt

REAL(8)::at0,gnt0,ctt0,ht0,pt0,wt0,bt0,taut0,tauvalt0,rt0,fmt0  
REAL(8),DIMENSION(nt2)::at2,gnt2,ctt2,ht2,cnt2,pt2,rert2,wt2
REAL(8),DIMENSION(nt2)::bt2,rt2,taut2,tauvalt2,yt2,ytv2,ynt2,ytt2,yttv2,ct2,fmt2
        
REAL(8),DIMENSION(nsam)::corgh,corgy,corgp,corgc,corgr,corga,corgb,corpgp,corgyNyT
REAL(8),DIMENSION(nsam)::corar,corha,coryh,corypr,coryc,coryb,corab
REAL(8),DIMENSION(nsam)::coryr,meanr,meangyn,meangy,corgtau,corgtauval,corgby,coryy,coryyn

REAL(8),DIMENSION(nsam)::meanyty,meanby,meanbyss,meanbyt,meanw,meanp,meanh
REAL(8),DIMENSION(nsam)::meanh1,meanb,meanyt2,meang,meanfm
  	
REAL(8),DIMENSION(nsam)::sdy,sdyr,sdyn,sdpnyn,sdyt,sdc,sdct,sdgyn,sdgy,sdg,sdg2,sdg3,sdp,sdr,sdcsdy,sdunemp,sdrer
REAL(8),DIMENSION(nsam)::sdctsdy,meantau,meantauval,meanrev,sdb,sdtau,sdtau2,sdtauval,sddtauval

REAL(8)::scorgh,scorgy,scorgp,scorgc,scorgr,scorga,scorgb,scorar,scorha,scorgby,scorgtau,scorgtauval,scoryb
REAL(8)::scorab,scoryy,scoryyn,scorpgp,scorgyNyT
REAL(8)::scoryh,scorypr,scoryc,scoryr,smeanr,smeangyn,smeangy,smeang,smeanb

REAL(8)::smeanyty,smeanyt2,smeanby,smeanbyt,smeanbyss,smeanw,smeanp,smeanh,smeanh1,smeantau,smeantauval,smeanrev,smeanfm  	
REAL(8)::ssdy,ssdyr,ssdyn,ssdyt,ssdc,ssdgyn,ssdgy,ssdp,ssdr,ssdg,ssdg2,ssdg3,ssdct,ssdcsdy,ssdctsdy,ssdtau,ssdtau2,ssdtauval,ssdb,ssdgtau,sfreq,ssdunemp,ssdpnyn,ssdrer,ssddtauval	

REAL(8)::auxpb(na2,nb2),dummy1(nt+n0),ones(na2,nb2),gnaut1aux(na2,nb2),minb

INTEGER::i_at(nt),i_bt(nt),i_at0(nt+n0),i_a0,bort0,acct0,acct02,deft0,i_at2(nt2)
INTEGER::i,j,ij,ii,jj,kk,jj3,k,t,i_a,i_b,ibpmin,ibpmax,gnmine,gnmaxe,i_aut0(1)
INTEGER::k1aux,k1,ksim,N_d,N_dd,aux,ia0,i0debt

REAL(8),DIMENSION(200)::auxa,auxgn,auxy,auxb

REAL(8),DIMENSION(nsam)::corgyr,corycr,coryrr,corgpr,corgbr
REAL(8)::scorgyr,scorycr,scoryrr,scorgpr,scorgbr

REAL(8)::deftemp(nt),bortemp(nt),acctemp(nt),q1c
REAL(8)::freqminb(nmc),freqmaxb(nmc),sfreqminb,sfreqmaxb  
REAL(8)::fm1(na2,nb2),fmn1(na2,nb2),freq(nmc)

REAL(8),DIMENSION(9,1000)::gnt_sim,ctt_sim,ht_sim,taut_sim,pt_sim,rt_sim,bt_sim,wt_sim,at_sim,deft_sim,bort_sim,acct_sim  
REAL(8),DIMENSION(9,1000)::gnt_simc,ctt_simc,ht_simc,taut_simc,pt_simc,rt_simc,bt_simc,wt_simc,at_simc  

INTEGER::bc_num,choice_fm,choice_sim2,choice_distfm,choice_data
REAL(8)::c_drop_def,h_drop_def,cn_def(2,1000),c_def(2,1000),h_def,h_bdef 
REAL(8)::freqbc(nsam),sfreqbc,tauval00

CHARACTER(len=26)::filename

1002 FORMAT(1000F7.2)
     
fm1=0d0
choice_fm=0         !compute MC stats for fiscal multipliers
choice_sim2=0       !compute default dynamics and cons drop at default
choice_distfm=0     !compute asymptotic distribution of fiscal multipliers

IF (choice_ols.EQ.1) THEN
    OPEN(55,FILE='output//ols//timeseries.txt')
    OPEN(56,FILE='output//ols//timeseriesL.txt')
    OPEN(57,FILE='output//ols//timeseriesH.txt')
ENDIF

IF (choice_fm.EQ.1) THEN
    OPEN(1,FILE='output//fm.txt')
            
    DO j = 1,nb2
        DO i = 1,na2
        READ(1,*) fm1(i,j)
    ENDDO
    ENDDO
 
CLOSE(1)

ENDIF

IF (export_data.EQ.1) THEN     
    OPEN(223,FILE='output//seriesb_are.txt')  
    OPEN(230,FILE='output//seriesia_are.txt')    
ENDIF

seed=1234658753
i_a0 = 1
minb = 0d0
k=1
ksim   = 0          !number of samples to generate around default episodes

IF (choice_dist .EQ. 1) THEN   
   OPEN(122,FILE='output//mcdistA2.txt',STATUS='REPLACE')    
ENDIF

choice_data=0

DO j=1,nmc
    CALL TIME_SERIES(n0+nt,na2,prob,DEXP(agrid2),DEXP(agrid2(i_a0)),i_at0)
    CALL r8vec_uniform_01(n0+nt,seed,dummy1)    !reentry
    
    i_b = LOCATE2(bgrid2,0d0) 
    
    acct0 = 1   !access in initial period
    deft0 = 0
    bort0 = 1
    tauval00=16d0
 
    DO i=1,n0+nt

        IF (i_at0(i).EQ.0) THEN
            PRINT*,'Problem here MC Stats'
            PRINT*,i,j,i_at0(i)
            PAUSE
            STOP
        ENDIF
        
        i_a = i_at0(i)
        acct02 = acct0
        deft0  = 0
        bort0  = 1

        IF ((acct0.EQ.0).AND.(dummy1(i)>phi)) THEN
            bort0 = 0 !bort(i)=0
            acct0 = 0 !acct(i+1)=0
        ELSEIF ((acct0.EQ.0).AND.(dummy1(i).LE.phi)) THEN
            bort0 = 1 !bort(i)=0
            acct0 = 1 !acct(i+1)=0
        ENDIF
        
        IF (acct0.EQ.1) THEN 
            IF (def1(i_a,i_b).EQ.1d0) THEN
                bort0 = 1 !bort(i)=1
                acct0 = 1 !acct(i+1)=1
            ELSE
                bort0 = 0  !bort(i)=0
                acct0 = 0  !acct(i+1)=0
                deft0 = 1  !deft(i)=1
            ENDIF
        ENDIF

        IF (bort0.EQ.1) THEN    !bort(i)
            at0  = DEXP(agrid2(i_a))
            gnt0 = gn1(i_a,i_b)
            ctt0 = ct1(i_a,i_b)
            ht0  = h1(i_a,i_b)
            pt0  = p1(i_a,i_b)
            wt0  = w1(i_a,i_b)
            taut0= tau1(i_a,i_b)
            tauvalt0= tauval1(i_a,i_b)
            fmt0 = fm1(i_a,i_b)

            i_b   = bp1(i_a,i_b)
            bt0   = bgrid2(i_b)/(lambda+r)
            acct0 = 1
            rt0   = (1d0+q1(i_a,i_b)**(-1d0)-lambda)/(1d0+r)-1d0       

        ELSEIF (bort0.EQ.0) THEN
            at0  = DEXP(agrid2(i_a))
            gnt0 = gnaut1(i_a)
            ctt0 = ctaut1(i_a)
            ht0  = haut1(i_a)
            pt0  = paut1(i_a)
            wt0  = waut1(i_a)
            taut0= tauaut1(i_a)
            tauvalt0= tauvalaut1(i_a)
            fmt0 = -1d0 
            rt0  = -1d0
            i_b  = i_bzero2
        ENDIF
    
        IF (i.GT.n0) THEN
            i_at(i-n0)= i_a

            at(i-n0)  = at0
            gnt(i-n0) = gnt0
            ctt(i-n0) = ctt0
            ht(i-n0)  = ht0
            pt(i-n0)  = pt0
            wt(i-n0)  = wt0
            i_bt(i-n0)= i_b
            rt(i-n0)  = rt0
            taut(i-n0)= taut0
            tauvalt(i-n0)= tauvalt0
            dtauvalt(i-n0)=(tauvalt0/tauval00-1d0)*100
            fmt(i-n0) = fmt0            
            deft(i-n0)= deft0
            acct(i-n0)= acct02
            bort(i-n0)= bort0
        ENDIF
        
        tauval00=tauvalt0       !previous-period tax rate
    ENDDO  !i

    bt   = bgrid2(i_bt)/(lambda+r)
    minb = MIN(minb,MINVAL(bt) )
    
    IF (mu .NE. 0d0) THEN
        rert = (omegac**(1d0/(1d0+mu))+ (1d0-omegac)**(1d0/(1d0+mu))*pt**(mu/(1d0+mu)))**(-(1d0+mu)/mu)
    ELSE
        rert = pt**(-(1d0-omegac))
    ENDIF

    IF (choice_dist.EQ.1) THEN
        i=nt
        WRITE(122,'(I7)') na2*nb2*(1-INT(bort(i)))+INT(bort(i))*(i_bt(i)-1)*na2 + i_at(i)
    ENDIF
    
    IF ((j.EQ.10).OR.(choice_data.EQ.1)) THEN
        OPEN(1,FILE='output//seriesp.txt')
        OPEN(2,FILE='output//seriesh.txt')
        OPEN(3,FILE='output//seriesb.txt')
        OPEN(4,FILE='output//seriesct.txt')
        OPEN(5,FILE='output//seriesr.txt')
        OPEN(7,FILE='output//seriesgn.txt')
        OPEN(8,FILE='output//seriestau.txt')    
        OPEN(88,FILE='output//seriestauval.txt')    
        OPEN(9,FILE='output//seriesa.txt')    
        OPEN(11,FILE='output//seriesbor.txt')    
        DO i=1,(1-choice_data)*205+choice_data*16
    
            WRITE(1,*) pt(i)
            WRITE(2,*) ht(i)
            WRITE(3,*) bt(i)
            WRITE(4,*) ctt(i)
            WRITE(5,*) rt(i)*100d0
            WRITE(7,*) gnt(i)
            WRITE(8,*) taut(i)
            WRITE(88,*) tauvalt(i)
            WRITE(9,*) at(i)
            WRITE(11,*) bort(i)
        ENDDO
    
        CLOSE(1)
        CLOSE(2)
        CLOSE(3)
        CLOSE(4)
        CLOSE(5)
        CLOSE(7)
        CLOSE(8)
        CLOSE(88)
        CLOSE(9)
        CLOSE(11)
    
    PRINT*, ' '
    PRINT*, 'TIME SERIES EXPORTED'
    PRINT*, ' '
    
    IF (choice_data.EQ.1) RETURN
    ENDIF

    IF (DBLE(SUM(acct)).GT.0d0) freq(j)= 100d0*(1d0-(1d0-DBLE(SUM(deft))/DBLE(SUM(acct))))
    
    IF (psi1.GE.500d0) THEN
        freqminb(j) = 100d0*DBLE(COUNT(bgrid2(i_bt).LE.bmin+1d-3))/DBLE(nt)
        freqmaxb(j) = 100d0*DBLE(COUNT(bgrid2(i_bt).GE.bmax-1d-3))/DBLE(nt)
    ENDIF
    
    IF (j.EQ.1) PRINT*,j,'Default frequency = ',freq(j),SUM(deft),SUM(acct)

!*******************************************************
!Compute dynamics around default events
    N_d = 4         !# periods of autarky after default event
    N_dd= 9         !# periods of window around default event
    
    IF (choice_sim2.EQ.1) THEN      !compute simulation paths

    ii=Nt2

    DO WHILE (ksim .LT. 1000)     !1000 simulations around default events

	    deftemp = deft(:)
	    bortemp = bort(:)
        acctemp = acct(:)
	    aux = 0                         !#subsamples of default events in the simulation
    
        DO jj=0,nt-(N_dd+2)            !index for time t

        ii = nt-jj-aux*(N_dd+2)-N_d   !2 quarters before w/access, i: index for default event 

        IF (ii .LE. (N_d+2+1)) THEN 
            GO TO 20202
        ELSEIF (ii .GT. (N_d+2+1)) THEN

	        IF ((bortemp(ii-1).EQ.1) .AND. (bortemp(ii).EQ.0)) THEN !default at time i

		    IF ((MINVAL(acctemp(ii-(N_d+2):ii-1)).EQ.1).AND.(MAXVAL(bortemp(ii+1:ii+N_d)).EQ.0)) THEN !no default in previous N_d+2 periods, autarky in following N_d periods

                ksim = ksim+1 !another sample
	    	    aux = aux+1   !another sample in this MC simulation

        	    IF (ksim.GT.1000) THEN
                    GO TO 20202 
    		    ENDIF

                WRITE(44,'(I7)') ii

                gnt_sim(:,ksim) = gnt(ii-N_d:ii+N_d)  
			    ctt_sim(:,ksim) = ctt(ii-N_d:ii+N_d)  
    		    ht_sim(:,ksim)  = ht(ii-N_d:ii+N_d)  
    		    taut_sim(:,ksim)= taut(ii-N_d:ii+N_d)  
			    pt_sim(:,ksim)  = pt(ii-N_d:ii+N_d)  
                rt_sim(:,ksim)	= rt(ii-N_d:ii+N_d)
    		    bt_sim(:,ksim)	= bt(ii-N_d:ii+N_d+1)
                wt_sim(:,ksim)	= wt(ii-N_d:ii+N_d)
    		    at_sim(:,ksim)  = at(ii-N_d:ii+N_d)  
    		    deft_sim(:,ksim)= deft(ii-N_d:ii+N_d+1)
                acct_sim(:,ksim)= acct(ii-N_d:ii+N_d)
    		    bort_sim(:,ksim)= bort(ii-N_d:ii+N_d)  
   
       		    at_simc(1:4,ksim)  = at(ii-N_d:ii-1) 
       		    ctt_simc(1:4,ksim) = ctt(ii-N_d:ii-1) 
       		    gnt_simc(1:4,ksim) = gnt(ii-N_d:ii-1) 
       		    ht_simc(1:4,ksim)  = ht(ii-N_d:ii-1) 
      		    taut_simc(1:4,ksim)= taut(ii-N_d:ii-1) 
                pt_simc(1:4,ksim)  = pt(ii-N_d:ii-1) 
                rt_simc(1:4,ksim)  = rt(ii-N_d:ii-1) 
                wt_simc(1:4,ksim)  = wt(ii-N_d:ii-1) 
                bt_simc(1:4,ksim)  = bt(ii-N_d:ii-1)
                
                i0debt = i_bt(ii-1)
            
                DO jj3=0,N_d    
   
                ia0 = i_at(ii+jj3)           
                at_simc(jj3+N_d+1,ksim)  = agrid(ia0)   
   
			    ctt_simc(jj3+N_d+1,ksim) = ct1(ia0,i0debt)  
			    gnt_simc(jj3+N_d+1,ksim) = gn1(ia0,i0debt)  
			    ht_simc(jj3+N_d+1,ksim)  = h1(ia0,i0debt)  
			    taut_simc(jj3+N_d+1,ksim)= tau1(ia0,i0debt)  
			    pt_simc(jj3+N_d+1,ksim)  = p1(ia0,i0debt)  
			    wt_simc(jj3+N_d+1,ksim)  = w1(ia0,i0debt)  
            
			    i0debt = bp1(ia0,i0debt)
			    bt_simc(jj3+N_d+1,ksim)  = bgrid2(i0debt)/(lambda+r)
            
                q1c = q1(ia0,i0debt)
                rt_simc(jj3+N_d+1,ksim) = (1d0+q1c**(-1d0)-lambda)/(1d0+r)-1d0
			             
                ENDDO
                
            ENDIF
	        ENDIF
        ENDIF
    
    ENDDO !jj

    IF ((j.EQ.Nmc).AND.(ksim.LT.1000)) THEN
	    PRINT*, 'Problem with Sampling: Insufficient series to generate 1000 samples'
	    PRINT*, 'Number of samples generated =',ksim
	    EXIT
    ENDIF

    ENDDO

    ENDIF !choice_sim2
!***********************************************
20202 t=1

DO WHILE ((t.LE.nt).AND.(k.LE.nsam))

IF ((t.GT.4).AND.(t+nt2-1.LE.nt)) THEN
    
    IF (MINVAL(bort(t-4:t+nt2-1)).EQ.1) THEN     !subsample qualifies: no default in next nt2 periods

        i_at2(:)= i_at(t:t+nt2-1)        
        at2(:)  = at(t:t+nt2-1)
        gnt2(:) = gnt(t:t+nt2-1)
        ctt2(:) = ctt(t:t+nt2-1)
        ht2(:)  = ht(t:t+nt2-1)
        pt2(:)  = pt(t:t+nt2-1)
        rert2(:)= rert(t:t+nt2-1)
        wt2(:)  = wt(t:t+nt2-1)
        bt2(:)  = bt(t:t+nt2-1)
        rt2(:)  = rt(t:t+nt2-1)
        taut2(:)= taut(t:t+nt2-1)
        tauvalt2(:)= tauvalt(t:t+nt2-1)

        IF (MAXVAL(rt2).GE.1d3) THEN
            PRINT*,'j',j
            PRINT*,' '
            PAUSE
            DO jj3=1,32
            PRINT*,jj3,rt2(jj3),bort(jj3)
            PAUSE
            ENDDO
        ENDIF
        
        cnt2 = ht2**alpha-gnt2
        
        !check frequency govt b.c. binds
        bc_num=0
        DO kk=1,nt2
            IF (DABS(pt2(kk)*gnt2(kk)-taut2(kk)-(ctt2(kk)-at2(kk))).LE.1d-5) THEN
                bc_num=bc_num+1
            ENDIF
        ENDDO

        IF (export_data.EQ.1) THEN
            WRITE(223,*) bt(t-1)*(lambda+r)
            DO ij=1,nt2
                WRITE(230,'(I7)') i_at2(ij)
            ENDDO
        ENDIF
        
IF (mu .NE. 0d0) THEN
    ct2 = (omegac*ctt2**(-mu)+ (1d0-omegac)*cnt2**(-mu))**(-1d0/mu)
ELSE
    ct2 = (ctt2**(omegac))*(cnt2**(1d0-omegac))
ENDIF

ynt2 = ht2**alpha
ytt2 = at2
yt2  = ytt2 + pt2*ynt2 


IF ((choice_ols.EQ.1).AND.(k.LE.200)) THEN
    DO i=2,nt2
        WRITE(55,'(F15.10,F15.10,F15.10,F15.10)') gnt2(i),ytt2(i),ynt2(i),bt2(i-1)
        IF (DLOG(ytt2(i)).LT.-0.27982d0 -1.21*bt2(i-1)*(lambda+r)) THEN
            WRITE(56,'(F15.10,F15.10,F15.10,F15.10)') gnt2(i),ytt2(i),ynt2(i),bt2(i-1)*(lambda+r)
        ELSE
            WRITE(57,'(F15.10,F15.10,F15.10,F15.10)') gnt2(i),ytt2(i),ynt2(i),bt2(i-1)*(lambda+r)
        ENDIF            
    ENDDO
    IF (k.EQ.200) THEN
        CLOSE(55)
        CLOSE(56)
        CLOSE(57)
    ENDIF
ENDIF

meanr(k)   = mean(rt2)
meangyn(k) = mean(gnt2/ynt2)
meang(k)   = mean(gnt2)
meangy(k)  = mean(pt2*gnt2/yt2)
meangyn(k) = mean(gnt2/ynt2)

meanyty(k) = mean(ytt2/yt2)
meanyt2(k) = mean(yt2)

meanby(k)  = mean(bt2/yt2)

meanbyt(k) = mean(bt2/ytt2)
meanw(k)   = mean(wt2)
meanb(k)   = mean(bt2)
meanp(k)   = mean(pt2)
meanh(k)   = mean(ht2)

meanh1(k)  = SUM(-DBLE((ht2.EQ.1d0)))/DBLE(nt2)

freqbc(k)  = 100d0*DBLE(bc_num)/DBLE(nt2) 
meanfm(k)  = mean(fmt2)

meantau(k) = mean(taut2/yt2)
meantauval(k) = mean(tauvalt2)
meanrev(k) = mean(taut2)

sdyt(k)    = std(DLOG(ytt2))
sdy(k)     = std(DLOG(yt2))
sdyn(k)    = std(DLOG(ynt2))
sdpnyn(k)  = std(pt2*ynt2)
sddtauval(k)= stdabs(dtauvalt)

sdyr(k)    = std(DLOG(ytt2+meanp(k)*ynt2))
sdtau(k)   = std(DLOG(taut2))/std(DLOG(ytt2+meanp(k)*ynt2))
sdtau2(k)  = std(DLOG(taut2))/std(DLOG(yt2))
sdtauval(k)= std(tauvalt2)
sdc(k)     = std(DLOG(ct2))
sdunemp(k) = std(1d0-ht2)

sdgyn(k)   = std(gnt2/ynt2)
sdgy(k)    = 100d0*std(gnt2/yt2)
sdg(k)     = std(DLOG(meanp(k)*gnt2)) /std(DLOG(ytt2+meanp(k)*ynt2))      
sdg2(k)    = std(gnt2/ynt2)       
sdg3(k)    = std(DLOG(gnt2))                        
    
sdp(k)     = std(DLOG(pt2))
sdrer(k)   = std(DLOG(rert2))
sdr(k)     = 100d0*std(rt2)
sdcsdy(k)  = std(DLOG(ct2))/std(DLOG(ytt2+meanp(k)*ynt2))  
sdb(k)     = std(bt2)
sdct(k)    = std(DLOG(ct2))

sdctsdy(k) = std(DLOG(ctt2))/std(DLOG(yt2))  	

corga(k) = corr(DLOG(gnt2),DLOG(at2))
corgh(k) = corr(DLOG(gnt2),1d0-ht2)
corgy(k) = corr(DLOG(gnt2),DLOG(yt2))
corgp(k) = corr(DLOG(gnt2),DLOG(rert2))
corpgp(k)= corr(DLOG(pt2*gnt2),DLOG(rert2))
corgc(k) = corr(DLOG(gnt2),DLOG(ct2))
corgb(k) = corr(DLOG(gnt2),bt2)
corgr(k) = corr(DLOG(gnt2),rt2)

corgtau(k) = corr(DLOG(gnt2),DLOG(taut2))
corgtauval(k) = corr(DLOG(gnt2),tauvalt2)
corgby(k)  = corr(DLOG(gnt2),bt2/yt2)

corha(k) = corr(DLOG(at2),1d0-ht2)
corar(k) = corr(DLOG(at2),rt2)

coryyn(k)= corr(DLOG(ytt2),DLOG(pt2*ynt2))
coryy(k) = corr(DLOG(ytt2),DLOG(ynt2))
coryh(k) = corr(DLOG(ytt2+meanp(k)*ynt2),1d0-ht2)

coryc(k) = corr(DLOG(yt2),DLOG(ct2))
coryr(k) = corr(DLOG(yt2),rt2)
coryb(k) = corr(DLOG(yt2),bt2)
corab(k) = corr(DLOG(at2),bt2)

corgyr(k) = corr(DLOG(meanp(k)*gnt2),DLOG(ytt2+meanp(k)*ynt2))
corycr(k) = corr(DLOG(ytt2+meanp(k)*ynt2),DLOG(ct2))
coryrr(k) = corr(DLOG(ytt2+meanp(k)*ynt2),rt2)
corypr(k) = corr(DLOG(ytt2+meanp(k)*ynt2),DLOG(rert2))
corgpr(k) = corr(DLOG(meanp(k)*gnt2),DLOG(rert2))
corgbr(k) = corr(DLOG(meanp(k)*gnt2),bt2)
corgyNyT(k) = corr(gnt2/ynt2,DLOG(ytt2))

k = k+1
t = t+nt2-1+4 - 2000*(psi1.GT.500d0)

ELSEIF (MINVAL(bort(t-4:t+nt2-1)).EQ.0) THEN

    i_aut0 = MINLOC(bort(t-4:t+nt2-1))
    i_aut0 = i_aut0(1)+t       !time index for default episode         
    i_aut0 = i_aut0(1)+4       !4 periods after default

    t = i_aut0(1)
ENDIF

ELSE
    t=t+1
ENDIF

ENDDO  !WHILE        

ENDDO

IF (k.LT.nsam) THEN 
    PRINT*,'Insufficient number of access samples, k =',k
    PAUSE
ENDIF
    
IF (choice_dist.EQ.1) THEN    
    CLOSE(122) 
    PRINT*, ' '
    PRINT*, 'ERGODIC DISTRIBUTION SAVED'
    PRINT*, ' '
ENDIF


if (export_data.EQ.1) THEN
    CLOSE(223)
    CLOSE(230)
ENDIF


scorga = mean(corga) 
scorgh = mean(corgh) 
scorgy = mean(corgy) 
scorgyr= mean(corgyr) 

scorgb = mean(corgb) 
scorgp = mean(corgp) 
scorpgp= mean(corpgp) 
scorgc = mean(corgc)
scorgr = mean(corgr)
scorgtau = mean(corgtau)
scorgtauval = mean(corgtauval)
scorgby  = mean(corgby)
scorgyNyT = mean(corgyNyT)

scorha = mean(corha)
scorar = mean(corar) 

scoryh = mean(coryh)
scorypr= mean(corypr) 
scoryc = mean(coryc) 
scoryr = mean(coryr) 

scoryy = mean(coryy)
scoryyn= mean(coryyn)
scorycr= mean(corycr)
scoryrr= mean(coryrr)
scorgpr= mean(corgpr)
scorgbr= mean(corgbr) 
scoryb = mean(coryb) 
scorab = mean(corab) 

sfreq    = mean(freq) 
sfreqbc  = mean(freqbc)

sfreqminb= mean(freqminb)
sfreqmaxb= mean(freqmaxb) 
        
smeanr   = 100d0*mean(meanr)
smeang   = 100d0*mean(meang) 
smeangy  = 100d0*mean(meangy) 
smeangyn = 100d0*mean(meangyn) 

smeangyn = 100d0*mean(meangyn)

smeanyty = 100d0*mean(meanyty)
smeanby  = 100d0*mean(meanby)
smeanbyt = 100d0*mean(meanbyt)
smeanbyss= 100d0*mean(meanb/meanyt2)
smeanyt2 = 100d0*mean(meanyt2)

smeanyty = 100d0*mean(meanyty)
smeanby  = 100d0*mean(meanby)

smeanb   = mean(meanb) 
smeanw   = mean(meanw)
smeanp   = mean(meanp)
smeanh   = mean(meanh)
smeanh1  = mean(meanh1)
smeantau = mean(meantau)
smeantauval = mean(meantauval)
smeanrev = mean(meanrev)
smeanfm  = mean(meanfm)

ssdy     = mean(sdy)
ssdyt    = mean(sdyt)
ssdyn    = mean(sdyn)
ssdpnyn  = mean(sdpnyn)
ssdyr    = mean(sdyr)
ssdc     = mean(sdc)
ssdgyn   = mean(sdgyn)
ssdgy    = mean(sdgy)
ssdg     = mean(sdg)
ssdg2    = mean(sdg2)
ssdg3    = mean(sdg3)
ssdp     = mean(sdp)
ssdrer   = mean(sdrer)
ssdr     = mean(sdr)
ssdunemp = mean(sdunemp)
ssdcsdy  = mean(sdcsdy)
ssdb     = mean(sdb)
ssdtauval= 100d0*mean(sdtauval)
ssddtauval= mean(sddtauval)
ssdtau   = mean(sdtau)
ssdtau2  = mean(sdtau2)
ssdgtau  = mean(sdg/sdtau)
ssdct    = mean(sdct)

ssdctsdy = mean(sdctsdy)  	

ones   = 1d0

ibpmax = INT(MAXVAL(def1*bp1))
ibpmin = INT(MINVAL(def1*bp1+(ones-def1)*nb))

meanb    = smeanb
h_drop_def = 0d0
h_def  = 0d0
h_bdef = 0d0

IF (choice_sim2.EQ.1) THEN
h_drop_def = mean(ht_sim(4,:)-ht_sim(1,:)) 
h_def  = mean(ht_sim(5,:))
h_bdef = mean(ht_sim(3,:))


OPEN(1,FILE='output//def_win//sim//median_gnt_sim.txt')
OPEN(2,FILE='output//def_win//sim//median_ctt_sim.txt')
OPEN(3,FILE='output//def_win//sim//median_ht_sim.txt')
OPEN(4,FILE='output//def_win//sim//median_pt_sim.txt')
OPEN(5,FILE='output//def_win//sim//median_bt_sim.txt')
OPEN(66,FILE='output//def_win//sim//median_rt_sim.txt')
OPEN(7,FILE='output//def_win//sim//median_acct_sim.txt')
OPEN(8,FILE='output//def_win//sim//median_yt_sim.txt')
        
DO j=1,9
    
    WRITE(1,*) quantile(gnt_sim(j,:),0.5d0)  
    WRITE(2,*) quantile(ctt_sim(j,:),0.5d0)  
    WRITE(3,*) quantile(ht_sim(j,:),0.5d0)  
    WRITE(4,*) quantile(pt_sim(j,:),0.5d0)  
    WRITE(5,*) quantile(bt_sim(j,:),0.5d0)  
    WRITE(66,*) quantile(rt_sim(j,:),0.5d0)  
    WRITE(7,*) quantile(bort_sim(j,:),0.5d0)  
    WRITE(8,*) quantile(at_sim(j,:),0.5d0)  
ENDDO

CLOSE(1)
CLOSE(2)
CLOSE(3)
CLOSE(4)
CLOSE(5)
CLOSE(66)
CLOSE(7)
CLOSE(8)

OPEN(51,FILE='output//def_win//sim//median_gnt_simc.txt')
OPEN(52,FILE='output//def_win//sim//median_ctt_simc.txt')
OPEN(53,FILE='output//def_win//sim//median_ht_simc.txt')
OPEN(54,FILE='output//def_win//sim//median_pt_simc.txt')
OPEN(55,FILE='output//def_win//sim//median_bt_simc.txt')
OPEN(56,FILE='output//def_win//sim//median_rt_simc.txt')
OPEN(57,FILE='output//def_win//sim//median_yt_simc.txt')

DO j=1,9
    WRITE(51,*) quantile(gnt_simc(j,:),0.5d0)  
    WRITE(52,*) quantile(ctt_simc(j,:),0.5d0)  
    WRITE(53,*) quantile(ht_simc(j,:),0.5d0)  
    WRITE(54,*) quantile(pt_simc(j,:),0.5d0)  
    WRITE(55,*) quantile(bt_simc(j,:),0.5d0)  
    WRITE(56,*) quantile(rt_simc(j,:),0.5d0)  
    WRITE(57,*) quantile(at_simc(j,:),0.5d0)  
ENDDO

CLOSE(51)
CLOSE(52)
CLOSE(53)
CLOSE(54)
CLOSE(55)
CLOSE(56)
CLOSE(57)

ENDIF

    
PRINT*, ' '
PRINT*, 'EXPORTING PARAMETERS AND MODEL''S STATISTICS TO OUTPUT.TXT'

WRITE(filename,'("output//final_",I3,".txt")') version_code
OPEN(1,FILE=filename,STATUS='REPLACE')

	WRITE(1,'(A40)') 'PARAMETERS:'
	WRITE(1,'(A40)') '=========='
	WRITE(1,'(2X)')     
	WRITE(1,*) 'PREFERENCES AND OUTPUT:'
	WRITE(1,'(A40,F10.6)') 'Discount factor, beta.............= ',beta*100d0
	WRITE(1,'(A40,F10.3)') 'Risk Aversion parameter, sigma....= ',sigma
	WRITE(1,'(A40,F10.3)') 'Risk Aversion parameter gN, sigmag= ',sigmag
	WRITE(1,'(A40,F10.3)') 'Tradable weight, omegac...........= ',omegac
	WRITE(1,'(A40,F10.3)') 'Elasticity of substitution, mu....= ',mu
	WRITE(1,'(A40,F10.3)') 'Labor share in nontradables, alpha= ',alpha
	WRITE(1,'(A40,F10.3)') 'Imported i share tradables, gamma.= ',gamma
	WRITE(1,'(A40,F10.5)') 'Utility parameter for gN, chi.....= ',chi
	WRITE(1,'(A40,F10.3)') 'Reentry probability, phi..........= ',phi
    WRITE(1,'(A40,F10.3)') 'Maturity parameter, lambda........= ',lambda	
	WRITE(1,'(A40,F10.3)') 'Iceberg cost parameter, icost.....= ',icost
    WRITE(1,'(A40,F10.3)') 'Constant tax rate, taxrate........= ',taxrate	
	WRITE(1,'(A40,F10.3)') 'Output cost parameter, yloss......= ',yloss	
    WRITE(1,'(A40,F10.5)') 'Output cost parameter, psi1.......= ',psi1	
    WRITE(1,'(A40,F10.5)') 'Output cost parameter, ps2........= ',psi2	
	WRITE(1,'(A40,F10.5)') 'Tax penalty parameter, taucost....= ',taucost
	WRITE(1,'(A40,I7)')    'Choice tax hat (=0 none)..........= ',choice_tax_hat
	WRITE(1,'(A40,F10.5)') 'Wage lower bound, wbar............= ',wbar
    WRITE(1,'(A40,F10.5)') 'Relative cons. cU/cE, alphac......= ',alphac
	WRITE(1,'(A40,F10.3)') 'Risk-free rate, r.................= ',r
	WRITE(1,'(2X)') 
	 
	WRITE(1,*) 'GRIDS:'
	WRITE(1,'(A40,I7)')    'Size grid tax.....................= ',ntau
	WRITE(1,'(A40,I7)')    'Size grid debt....................= ',nb
	WRITE(1,'(A40,I7)')    'Size grid TFP shock...............= ',na
	WRITE(1,'(A40,I7)')    'Number deviations TFP grid........= ',ndevs
    WRITE(1,'(A40,F10.3)') 'Minimum debt level, bmin..........= ',bmin
	WRITE(1,'(A40,F10.3)') 'Maximum debt level, bmax..........= ',bmax
	WRITE(1,'(A40,F10.3)') 'Lowest endog debt level...........= ',bgrid2(ibpmin)
	WRITE(1,'(A40,F10.3)') 'Highest endog debt level..........= ',bgrid2(ibpmax)
    WRITE(1,'(A40,F10.3)') 'Highest debt level simulation.....= ',minb
	WRITE(1,'(A40,I7)')    'Lowest endog debt index...........= ',ibpmin
	WRITE(1,'(A40,I7)')    'Highest endog debt index..........= ',ibpmax
    WRITE(1,'(A40,F10.2)') 'Fraction of debt grid in use......= ',DBLE((ibpmax-ibpmin+1)/nb)
    WRITE(1,'(A40,I7)')    'Monot V Rep (=1 really okay)......= ',monot_v0
    
IF (psi1.GE.-500d0) THEN
	WRITE(1,'(A40,F10.3)') 'freq(bmin hits).....................= ',sfreqminb 
	WRITE(1,'(A40,F10.3)') 'freq(bmax hits).....................= ',sfreqmaxb 
ENDIF
    
    WRITE(1,'(2X)') 
	WRITE(1,*) 'STOCHASTIC PROCESSES FOR PRODUCTIVITY:'
	WRITE(1,'(A40,F10.3)') 'mua...............................= ',mua
	WRITE(1,'(A40,F10.3)') 'rhoa..............................= ',rhoa
	WRITE(1,'(A40,F10.3)') 'sigmaa............................= ',sigmaa
	WRITE(1,'(2X)') 
	WRITE(1,*) 'VALUE FCN ITERATION AND MC SIMULATIONS:'
    WRITE(1,'(A40,F10.6)')  'Tolerance........................= ',tol
	WRITE(1,'(A40,I7)')     'Number of VFcn iters.............= ',niterv
	WRITE(1,'(A40,I7)')     'Number of VFcn iters used........= ',counter_global
	WRITE(1,'(A40,I7)')     'Number simulations, nmc..........= ',nmc
	WRITE(1,'(A40,I7)')     'Number access samples, nsam......= ',nsam
    WRITE(1,'(A40,I7)')     'Number periods, nt...............= ',nt
	WRITE(1,'(A40,I7)')     'No. periods access sample, nt2...= ',nt2
	WRITE(1,'(A40,I7)')     'Number burnt periods, nout.......= ',n0
	WRITE(1,'(A40,I7,A,I7,A)')  'Time to convergence..............= ',hours,'hs',mins,'mins'
	WRITE(1,'(2X)') 
	WRITE(1,'(2X)') 
	WRITE(1,'(2X)') 
	WRITE(1,'(2X)') 
	WRITE(1,'(A40)') 'MC STATISTICS:'
	WRITE(1,'(A40)') '=============='
	WRITE(1,'(2X)') 
	WRITE(1,'(A40,F10.3)') 'cor(gn,yT)........................= ',scorga 
	WRITE(1,'(A40,F10.3)') 'cor(gn,unemp).....................= ',scorgh  
	WRITE(1,'(A40,F10.3)') 'cor(gn,y).........................= ',scorgy 
	WRITE(1,'(A40,F10.3)') 'cor(gn,y real)....................= ',scorgyr 
    WRITE(1,'(A40,F10.3)') 'cor(gn,b).........................= ',scorgb 
	WRITE(1,'(A40,F10.3)') 'cor(gn real,b)....................= ',scorgbr 
	WRITE(1,'(A40,F10.3)') 'cor(gn,b/y).......................= ',scorgby 
	WRITE(1,'(A40,F10.3)') 'cor(gn,rer).......................= ',scorgp 
	WRITE(1,'(A40,F10.3)') 'cor(pngn,rer).....................= ',scorpgp 
    WRITE(1,'(A40,F10.3)') 'cor(gn,c).........................= ',scorgc  
	WRITE(1,'(A40,F10.3)') 'cor(gn,r).........................= ',scorgr 
	WRITE(1,'(A40,F10.3)') 'cor(gn,tau).......................= ',scorgtau  
	WRITE(1,'(A40,F10.3)') 'cor(gn,tax rate)..................= ',scorgtauval  
	WRITE(1,'(A40,F10.3)') 'cor(gn/yN,yT).....................= ',scorgyNyT  
	WRITE(1,'(2X)') 
	WRITE(1,'(A40,F10.3)') 'cor(yT,unemp).....................= ',scorha  
	WRITE(1,'(A40,F10.3)') 'cor(yT,yN)........................= ',scoryy 
	WRITE(1,'(A40,F10.3)') 'cor(yT,pN*yN).....................= ',scoryyn
	WRITE(1,'(A40,F10.3)') 'cor(y,h)..........................= ',scoryh   
	WRITE(1,'(A40,F10.3)') 'cor(y real,rer)...................= ',scorypr  
	WRITE(1,'(A40,F10.3)') 'cor(y,b)..........................= ',scoryb  
	WRITE(1,'(A40,F10.3)') 'cor(a,b)..........................= ',scorab  
    WRITE(1,'(A40,F10.3)') 'cor(y,c)..........................= ',scoryc  
    WRITE(1,'(A40,F10.3)') 'cor(y2 real,c)....................= ',scorycr   
	WRITE(1,'(A40,F10.3)') 'cor(y,r)..........................= ',scoryr  
	WRITE(1,'(A40,F10.3)') 'cor(y real,r).....................= ',scoryrr  
	WRITE(1,'(A40,F10.3)') 'cor(yT,r).........................= ',scorar  
	WRITE(1,'(2X)') 
	WRITE(1,'(A40,F10.3)') 'freq(default).....................= ',sfreq 
	WRITE(1,'(2X)') 
	WRITE(1,'(A40,F10.3)') 'mean(r)...........................= ',smeanr  	 	
	WRITE(1,'(A40,F10.3)') 'mean(y real)......................= ',smeanyt2  	
	WRITE(1,'(A40,F10.3)') 'mean(gn)..........................= ',smeang  
    WRITE(1,'(A40,F10.3)') 'mean(gn/y)........................= ',smeangy  
    WRITE(1,'(A40,F10.3)') 'mean(gn/yn).......................= ',smeangyN  	 	
	WRITE(1,'(A40,F10.3)') 'mean(yt/y)........................= ',smeanyty  	
	WRITE(1,'(A40,F10.3)') 'mean(b/y).........................= ',smeanby  	
	WRITE(1,'(A40,F10.3)') 'mean(b/yssmc).....................= ',smeanbyss  		
	WRITE(1,'(A40,F10.3)') 'mean(b)...........................= ',smeanb  	 	
	WRITE(1,'(A40,F10.3)') 'mean(w)...........................= ',smeanw  
	WRITE(1,'(A40,F10.3)') 'mean(rer).........................= ',smeanp  	
	WRITE(1,'(A40,F10.3)') 'mean(h)...........................= ',smeanh  	
	WRITE(1,'(A40,F10.3)') 'mean(tau/y).......................= ',smeantau  	
	WRITE(1,'(A40,F10.3)') 'mean(tax rate)....................= ',smeantauval  	
	WRITE(1,'(A40,F10.3)') 'mean(revenues)....................= ',smeanrev 
	WRITE(1,'(A40,F10.3)') 'mean(fiscal multi)................= ',smeanfm  	    
	WRITE(1,'(A40,F10.3)') 'mean(time with h=1)...............= ',smeanh1  	
	WRITE(1,'(A40,F10.3)') 'mean(h drop at default)...........= ',h_drop_def
	WRITE(1,'(A40,F10.3)') 'mean(h at default)................= ',h_def
	WRITE(1,'(A40,F10.3)') 'mean(h before default)............= ',h_bdef
	WRITE(1,'(A40,F10.3)') 'freq(govt b.c. binds).............= ',sfreqbc     
	WRITE(1,'(2X)') 	
	WRITE(1,'(A40,F10.3)') 'std.dev.(y real)..................= ',ssdyr  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(y).......................= ',ssdy 
	WRITE(1,'(A40,F10.3)') 'std.dev.(yt)......................= ',ssdyt 
    WRITE(1,'(A40,F10.3)') 'std.dev.(yn)......................= ',ssdyn 
    WRITE(1,'(A40,F10.3)') 'std.dev.(pn*yn2)..................= ',ssdpnyn 
	WRITE(1,'(A40,F10.3)') 'std.dev.(pg)/std.dev.(tau)........= ',ssdgtau  	
    WRITE(1,'(A40,F10.3)') 'std.dev.(pg)/std.dev.(y real).....= ',ssdg 
	WRITE(1,'(A40,F10.3)') 'std.dev.(c).......................= ',ssdc  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(tau)/std.dev.(y real)....= ',ssdtau  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(tax rate)................= ',ssdtauval  	
	WRITE(1,'(A40,F10.3)') 'abs.dev.(tax rate)................= ',ssddtauval  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(gn)......................= ',ssdg3  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(gn/yn)...................= ',ssdg2  		
    WRITE(1,'(A40,F10.3)') 'std.dev.(b).......................= ',ssdb 	
	WRITE(1,'(A40,F10.3)') 'std.dev.(p).......................= ',ssdp  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(rer).....................= ',ssdrer 
    WRITE(1,'(A40,F10.3)') 'std.dev.(r).......................= ',ssdr  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(1-h).....................= ',ssdunemp  	
	WRITE(1,'(A40,F10.3)') 'std.dev.(ct)......................= ',ssdct  	 
  
    CLOSE(1)
    
    
xvector(1:6)=0d0

xvector(1) =smeanby
xvector(2) =smeangy
xvector(3) =smeanr
xvector(4) =ssdr
xvector(5) =ssdgtau
xvector(6) =mean(meanh)

PRINT*, ' '
PRINT*, 'EXPORTING PARAMETERS AND MODEL''S STATISTICS TO OUTPUT.TXT'

RETURN
END SUBROUTINE MC_SIM


